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Abstract 

OpenProp  is  an  open  source  propeller  and  turbine  design  and  analysis  code  that  has  been  in 
development  since  2007  by  MIT  graduate  students  under  the  supervision  of  Professor  Richard 
Kimball.  In  order  to  test  the  performance  predictions  of  OpenProp  for  axial  flow  hydrokinetic 
turbines,  a  test  fixture  was  designed  and  constructed,  and  a  model  scale  turbine  was  tested.  Tests 
were  conducted  in  the  MIT  water  tunnel  for  tip  speed  ratios  ranging  from  1.55  to  7.73. 

Additional  code  was  also  written  and  added  to  OpenProp  in  order  to  implement  ABS  steel 
vessels  rules  for  propellers  and  calculate  blade  stress.  The  blade  stress  code  was  used  to  conduct 
a  fatigue  analysis  for  a  model  scale  propeller  using  a  quasi-steady  approach. 

Turbine  test  results  showed  that  OpenProp  provides  good  performance  predictions  for  the  on- 
design  operational  condition  but  that  further  work  is  needed  to  improve  performance  predictions 
for  the  off-design  operational  condition.  Fatigue  analysis  results  show  that  reasonable  estimates 
of  propeller  blade  fatigue  life  can  be  obtained  using  a  relatively  simple  method.  Calculated  blade 
stress  distributions  agree  with  previously  published  data  obtained  with  more  sophisticated  and 
time  consuming  calculation  techniques. 
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Introduction 


Since  2007,  graduate  students  at  the  Massachusetts  Institute  of  Technology  (MIT)  have  been 
developing  an  open  source  propeller  and  turbine  design  and  analysis  tool  under  the  supervision 
of  Professor  Richard  Kimball.  The  tool  is  a  set  of  open  source  MATLAB®  scripts  published 
under  the  GNU  General  Public  License  which  are  capable  of  performing  design  and  analysis 
studies  for  open  and  ducted  propellers  as  well  as  axial  flow  turbines.  This  suite  of  MATLAB® 
scripts  is  called  OpenProp.  OpenProp  propeller  design  capabilities  include  performing 
parametric  studies  of  propellers  using  various  propeller  diameters,  number  of  blades  and  rotation 
speeds.  Propeller  analysis  features  include  perfonning  off-design  and  cavitation  analyses.  A  gap 
in  OpenProp  capabilities  was  the  inability  to  evaluate  the  structural  adequacy  of  a  propeller  or 
turbine.  This  project  added  two  new  modules.  One  module  implements  American  Bureau  of 
Shipping  (ABS)  steel  vessel  rules  for  propellers  and  the  other  calculates  the  blade  surface  stress. 

Validation  of  OpenProp  turbine  and  propeller  performance  predictions  is  limited.  The  portion  of 
the  code  suite  which  designs  ducted  propellers  has  been  validated  against  the  US  Navy’s 
Propeller  Lifting  Line  (PLL)  code  with  excellent  correlation.  Several  experiments  have  been 
done  to  validate  open  propeller  perfonnance  predictions  using  a  modified  trolling  motor 
apparatus.  One  test  had  been  performed,  with  limited  success,  of  an  axial  flow  turbine.  No  tests 
had  been  performed  for  ducted  propellers.  Because  of  this  lack  of  experimental  validation  of 
OpenProp,  it  became  necessary  to  design  and  construct  a  propeller  and  turbine  test  fixture  that  is 
robust  and  can  easily  be  used  to  test  open  and  ducted  propellers  as  well  as  turbines.  This  project 
provided  a  test  fixture,  funded  by  MIT  SeaGrant,  which  can  be  used  in  a  water  tunnel  or  tow  tank 
to  provide  experimental  performance  results  which  can  be  used  to  validate  OpenProp 
performance  predictions. 

OpenProp  implements  the  vortex  lifting  line  method  to  quickly  achieve  a  propeller  or  axial  flow 
turbine  design.  The  lifting  line  method  of  propeller  design  has  some  limitations  but  is  an 
excellent  method  to  obtain  an  initial  design  which  can  be  refined  using  more  sophisticated  design 
techniques.  In  the  spirit  of  providing  initial  design  estimates,  this  project  also  completed  a  quasi¬ 
steady  fatigue  analysis  and  predicted  the  fatigue  life  of  a  propeller. 

This  paper  presents  the  results  of  testing,  blade  stress  calculations  and  fatigue  analysis. 
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Chapter  1  - Development ,  Capability  and  Limitations  of  OpenProp 


Development  of  OpenProp 

OpenProp  had  its  genesis  in  a  code  called  MATLAB®  Propeller  Vortex  Lattice  (MPVL)  which 
was  a  code  which  added  graphical  user  interfaces  to  Propeller  Vortex  Lattice  (PVL)  which  was 
developed  by  Kerwin  (2007)  for  his  propellers  course  at  the  Massachusetts  Institute  of 
Technology  (MIT).  Since  that  time  significant  capability  has  been  added  to  the  code  and 
additional  features  and  capability  are  being  developed. 

OpenProp  uses  a  lifting  line  method  to  model  blade  circulation  (Kerwin  2007).  The  lifting  line 
technique  has  been  well  established  and  was  implemented  by  Kerwin  for  preliminary  parametric 
propeller  design  for  the  US  Navy  in  a  code  called  PLL.  OpenProp  development  sought  to 
expand  and  enhance  the  capabilities  of  Kerwin’ s  code  and  make  the  software  more  user  friendly. 
A  full  explanation  for  the  theory  of  operation  of  OpenProp  has  been  given  by  Epps  (2010b). 


Capability  of  OpenProp 

A  table  showing  the  development  history  and  current  capability  of  OpenProp  is  shown  below. 


Date 

Event 

Persons 

Responsible 

Description 

2001 

PVL  Developed 

J.E.  Kerwin 

Lifting  line  design  code  used  for  Kerwin’ s 
propeller  class  at  MIT 

2007 

MPVL  Developed 
(Later  renamed 
OpenProp  vl.0) 

H.  Chung 

K.P.  D’Epagnier 

MATLAB®  version  of  PVL  which 
incorporated  GUIs  for  parametric  and  blade 
row  design  and  geometry  routines  for  CAD 
(Rhino)  interface.  This  code  used  a  Lerb’s 
criteria  optimizer.  Chung  (2007), 

D’Epagnier  (2007) 

2008 

Cavitation  Analysis 
Routines 

Developed 

C.J.  Peterson 

Using  Drela’s  XFOIL,  routines  and 
executables  were  developed  for  conducting 
propeller  cavitation  analysis,  Peterson 
(2008). 

2008 

OpenProp  v2.0 

J.M.  Stubblefield 

Added  capability  for  ducted  propeller  design 
Stubblefield  (2008) 

2009 

OpenProp  v2.3 

B.P.  Epps 

Added  new  optimizer.  Incorporated  routines 
of  Peterson,  added  off-design  analysis, 
corrected  errors  and  added  ability  to  design 
axial  flow  turbines  with  or  without  blade 
chord  optimization.  Theory  described  in 

Epps  (2010b). 

2010 

Contra-Rotating 
Propeller  Design 

D.  Laskos 

Added  the  capability  for  contra-rotating 
propeller  design  with  cavitation  analysis. 
Laskos (2010) 

Table  1:  Development  History  of  OpenProp 
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This  project  added  the  capability  to  calculate  blade  stress  and  implement  ABS  rules  for 
propellers.  Epps  continues  to  refine  and  expand  OpenProp  capabilities  and  is  currently  working 
on  codes  to  predict  propeller  perfonnance  during  shaft  reversals. 


Limitations  of  OpenProp 

OpenProp  uses  the  lifting  line  method  to  model  the  blade  circulation.  There  are  limits  in  regard 
to  using  this  method  in  propeller  design. 

1 .  Constant  Radius  Vortex  Helix  -  In  the  implementation  of  the  lifting  line  method,  the 
trailing  vorticity  is  assumed  to  be  of  constant  radius.  For  propellers,  it  is  known  that  the 
trailing  vorticity  helix  radius  actually  decreases.  This  simplification  has  been  made  to 
ease  the  complexity  of  calculating  the  influence  of  the  trailing  vorticity  on  the  blade  itself. 
The  errors  introduced  with  this  simplification  are  relatively  minor  as  shown  in  the 
experimental  data  comparison  in  this  paper  and  by  Stubblefield  (2008)  in  his  comparison 
of  OpenProp  predictions  to  a  more  established  propeller  design  code. 

2.  Skew  and  Rake  -  OpenProp  does  not  allow  the  designer  to  design  blades  with  skew  or 
rake 

3.  3D  Lifting  Surface  Effects  -  By  its  very  nature  a  design  tool  based  on  the  lifting  line 
method  cannot  account  for  the  3D  lifting  surface  effects.  The  result  is  that  the  translation 
from  calculated  hydrodynamic  characteristics  to  a  blade  geometry  which  produces  the 
desired  hydrodynamic  performance  is  difficult  and  in  most  cases  will  contain  errors. 
Therefore,  the  blade  geometry  that  is  generated  from  a  lifting  line  code  is  an 
approximation  and  one  should  not  expect  the  propeller  performance  measured  by 
experiment  to  completely  match  the  lifting  line  predicted  performance.  A  comparison 
and  discussion  of  predicted  perfonnance  versus  experimental  performance  for  U.S.  Navy 
Propeller  4119  is  contained  in  Epps  (2010a,  p.  224-225). 
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Chapter  2  -  Hydrokinetic  Turbine  Design  and  Construction 


In  propeller  design  the  overall  objective  is  to  generate  a  specified  thrust  while  minimizing  the 
torque  required  to  produce  it.  In  turbine  design  the  goal  is  to  maximize  the  torque  and  minimize 
the  thrust.  A  procedure  which  can  be  used  with  OpenProp  for  turbine  design  is: 

1.  Determine  expected  Cd  and  Cl.  Typical  ranges  for  these  quantities  are  0.008<Cd<0.03 
and  0.2<C/<0.5.  The  actual  values  for  these  parameters  are  dictated  by  the  choice  of 
blade  section  shape,  flow  regime  and  the  degree  of  blade  section  scaling. 

2.  Perform  parametric  design  study  using  expected  Cd/Cl  to  determine  number  of  blades 
and  tip  speed  ratio.  A  typical  value  for  this  ratio  is  0.06. 

3.  Select  a  design  point  from  the  parametric  study  of  step  2.  The  turbine  design  point  is 
characterized  by  the  number  of  turbine  blades  and  the  tip  speed  ratio.  In  general,  the 
more  blades  that  a  turbine  has  the  greater  its  efficiency.  However,  in  actual  application 
this  must  be  balanced  by  the  manufacture  costs  of  the  turbine. 

4.  Choose  the  turbine  diameter,  free  stream  flow  speed  and  rotation  rate  consistent  with  the 
chosen  tip  speed  ratio  in  step  3  above  such  that  desired  power  is  achieved.  Maximum 
turbine  diameter  is  dictated  by  the  water  depth  and  installation  scheme  where  the  turbine 
will  operate.  It  is  generally  desirable  to  maximize  the  turbine  diameter  in  order  to 
maximize  the  turbine’s  power  capacity.  Free  stream  flow  speed  is  determined  by  the 
flow  where  the  turbine  will  be  installed.  Desired  rotation  rate  will  be  effected  by  the 
electrical  generator  selected  for  use  with  the  turbine. 

5.  Perform  an  off-design  performance  analysis.  An  off-design  perfonnance  analysis  is 
necessary  to  obtain  an  overall  picture  of  the  time  average  power  that  the  turbine  will 
produce.  This  analysis  is  especially  important  for  tidal  turbines  where  there  is  a 
fluctuation  of  flow  speed. 

6.  Determine  the  span-wise  blade  chord  and  thickness  distribution.  This  step  is  where  the 
blade  geometry  is  detennined  to  produce  the  characteristics  determined  in  the  previous 
steps.  OpenProp  can  perform  this  step  automatically  by  using  the  chord  optimizer. 

7.  Perfonn  blade  stress  analysis.  A  blade  stress  analysis  is  necessary  to  ensure  the  structural 
adequacy  of  the  blades. 

The  above  procedure  was  used  to  design  the  turbine  which  was  tested  by  Epps.  Results  are 
shown  in  Epps  (2010a).  The  same  turbine  was  retested  as  part  of  this  project.  Step  7  was  not 
performed  as  part  of  the  design  process  because  the  stress  module  of  OpenProp  was  not  available 
at  that  time.  The  turbine  diameter  was  selected  as  the  maximum  diameter  which  could  be 
manufactured  using  the  available  rapid  prototyping  equipment  and  tested  in  the  water  tunnel  test 
section.  The  number  of  blades  was  also  dictated  by  the  desire  to  maximize  the  turbine  diameter 
and  two  blades  were  selected. 

Epps  (2010a)  describes  the  procedure  implemented  in  OpenProp  to  conduct  a  parametric  design 
study,  optimize  a  single  turbine  design  and  perform  off-design  analysis.  For  turbines  the  method 
can  be  summarized  as  setting  the  blade  circulation  less  than  zero  and  then  simultaneously  solving 
a  set  of  equations  such  that  the  resulting  variables  represent  a  physically  realizable  condition. 
Epps  (2009,  2010a)  also  discusses  the  correct  way  to  optimize  a  turbine  design. 
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Once  the  turbine  was  designed,  the  geometry  module  of  OpenProp  was  used  to  create  the  set  of 
points  that  represent  the  blade  surface  in  three  dimensions.  This  set  of  points  was  loaded  into 
SolidWorks®  via  a  macro  developed  for  this  purpose.  In  SolidWorks®,  the  blade  geometry  was 
turned  into  a  solid  which  was  multiplied  into  two  blades  and  attached  to  a  hub.  This  file  was 
saved  in  .stl  format  and  loaded  into  the  rapid  prototyping  machine  for  production.  The  model 
scale  turbine  that  was  generated  in  this  way  was  tested  by  Epps  and  as  part  of  this  project. 
Turbine  test  results  are  presented  in  the  next  section. 

Table  9.1  of  Epps  (2010a,  p.  280)  contains  the  turbine  design  parameters.  That  table  is 
reproduced  here  in  Table  2. 


Parameter 

Value 

Description 

Z 

2 

Number  of  blades 

n 

19.1  rev/s 

Rotation  rate 

D 

0.25  m 

Diameter 

Vs 

3  m/s 

Free  stream  speed 

Dhub 

0.08382 

Hub  diameter 

M 

20 

Number  of  panels 

P 

1000  kg/m3 

Water  density 

1 

5 

Tip  speed  ratio 

Climax 

0.5 

Maximum  allowable  lift  coefficient 

Table  2:  Key  Turbine  Parameters 


Figure  1:  Turbine  Drawing 
Courtesy  of  Epps 


Figure  2:  Turbine  Test 
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Chapter  3  -  Test  Procedure ,  Results  and  Comparison 
Test  Procedure 

Turbine  testing  was  performed  on  a  test  fixture  specifically  designed  for  this  purpose.  The  test 
fixture  was  funded  via  MIT  SeaGrant  and  its  design  and  construction  are  described  in  Chapter  7. 
Generally  the  test  procedure  consisted  of  measuring  the  shaft  torque  created  by  the  turbine  for 
various  rotation  rates  and  flow  speeds.  A  detailed  test  procedure  follows  with  test  points  in 
Table  3  and  test  results  in  Figure  9. 

Calibration 

Calibration  of  the  test  fixture  was  performed  by  hanging  known  weights  from  the  output  shaft  of 
the  test  fixture  and  reading  strain  gage  amplifier  output  voltage  using  Lab  View®.  This 
calibration  technique  is  a  static  calibration;  a  better  calibration  technique  for  this  type  of  testing 
would  have  been  a  dynamic  calibration.  However,  a  dynamic  calibration  is  more  complicated 
and  requires  additional  equipment  which  was  unavailable.  Lab  View®  was  connected  to  the  test 
fixture  in  an  identical  way  for  both  calibration  and  testing.  Results  of  the  calibration  are  shown 
Figure  3  and  Figure  4. 


LabView  Sensor  Thrust  Calibration 
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Figure  3:  Thrust  Calibration 
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Figure  4:  Torque  Calibration 


Because  the  motor  drive  used  for  these  tests  uses  pulse  width  modulation  (PWM)  at  300VDC 
and  because  the  signal  wires  are  running  alongside  the  power  cable  (inside  the  same  1.5  inch 
diameter  standpipe)  there  was  a  concern  that  the  Signal  to  Noise  Ratio  (SNR)  would  be  too  small 
to  be  able  to  effectively  measure  the  signal  voltage.  This  concern  was  allayed  by  performing 
spectral  analyses  on  the  measured  signal.  A  typical  result  of  these  analyses  is  shown  in  Figure  5. 
The  graph  shows  that  there  is  minimal  interference. 


Figure  5:  Spectrum  Analysis-600RPM,  1.69m/s 


Since  the  calibration  that  was  used  was  a  static  calibration,  it  is  necessary  to  correct  the  measured 
torque  with  the  friction  torque  in  order  to  determine  the  actual  torque  produced  by  the  turbine.  A 
graph  of  friction  torque  measured  at  various  rotation  rates  without  a  hub  or  turbine  attached,  but 
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with  the  test  fixture  submerged  in  the  test  section,  is  given  in  Figure  6.  These  values  are  used  to 
correct  the  torque  measured  by  the  sensor. 
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Figure  6:  Friction  Torque 
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Figure  7:  Measured  Torque  -  No  Correction 


Figure  7  shows  the  torque  measured  by  the  sensor  without  torque  correction.  The  data  points  of 
this  figure  represent  the  uncorrected  torque  values  which  were  measured  at  the  tip  speed  ratios 
listed  in  Table  3.  Comparison  of  Figure  6  and  Figure  7  shows  that  the  friction  torque  is  a 
relatively  small  value  compared  to  the  total  measured  torque. 


Test  Steps 

Testing  began  by  determining  the  tip  speed  ratios  that  would  bracket  the  turbine’s  design  point 
and  reproduce  the  entire  off-design  performance  curve  generated  by  OpenProp.  The  tip  speed 
ratios  which  were  used  in  this  test  are  shown  in  Table  3. 
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Tip  Speed  Ratio  -  X 


RPM 


200 

250 

300 

350 

400 

450 

500 

550 

600 

r%7i 

1.100 

2.38 

2.97 

3.57 

4.16 

4.76 

5.35 

5.95 

6.54 

7.14 

7.73 

1.185 

2.21 

2.76 

3.31 

3.87 

4.42 

4.97 

5.52 

6.08 

6.63 

7.18 

1269 

2.06 

2.58 

3.09 

3.61 

4.13 

4.64 

5.16 

5.67 

6.19 

6.70 

1.354 

1.93 

2.42 

2.90 

3.38 

3.87 

4.35 

4.83 

5.32 

5.80 

6.28 

1.438 

1.82 

2.28 

2.73 

3.19 

3.64 

4.10 

4.55 

5.01 

5.46 

5.92 

1.523 

1.72 

2.15 

2.58 

3.01 

3.44 

3.87 

4.30 

4.73 

5.16 

5.59 

1.608 

1.63 

2.04 

2.44 

2.85 

3.26 

3.66 

4.07 

4.48 

4.89 

5.29 

1.692 

1.55 

1.93 

2.32 

2.71 

3.09 

3.48 

3.87 

4.25 

4.64 

5.03 

Table  3:  Test  Tip  Speed  Ratios 


The  steps  taken  to  gather  the  data  displayed  in  Figure  9  are  outlined  below: 

1 .  Generate  Table  3  which  represents  the  test  points  at  which  data  was  gathered.  Flow 
speeds  selected  correspond  to  integer  speed  reference  number  increments  of  Figure  8 

2.  Set  water  tunnel  impeller  speed  to  create  desired  flow  speed  in  test  section 

3.  Command  desired  test  fixture  motor  rotation 

4.  Collect  torque  voltage  measurements  via  the  LabView®  interface.  Sample  rate  was 
set  at  500Hz.  Sample  time  was  5-10  seconds. 

5.  Increase  test  fixture  motor  rotation  rate 

6.  Wait  approximately  10  seconds  for  transient  behavior  to  subside 

7.  Repeat  steps  4-6  until  data  has  been  collected  for  every  rotation  rate  at  the  test  section 
flow  speed 

8.  Increase  test  section  flow  speed 

9.  Wait  approximately  one  minute  for  transient  behavior  to  subside. 

10.  Repeat  steps  3-9  until  all  data  has  been  collected. 


Conducting  the  test  in  the  order  listed  above  minimizes  the  time  required  to  collect  data  since  the 
transient  is  much  longer  for  a  water  tunnel  impeller  speed  change  than  for  a  test  fixture  motor 
speed  change. 

In  step  2,  the  water  flow  speed  in  the  tunnel  was  not  measured  directly.  Normal  mode  of 
operation  is  to  measure  the  flow  speed  in  the  test  section  using  a  Laser  Doppler  Velocimetry 
(LDV)  system;  however  the  LDV  system  was  not  operational  at  the  time  of  the  test.  Previous 
experimentation  in  the  water  tunnel  generated  Figure  8.  Figure  8  relates  impeller  rotation  rate  to 
test  section  flow  speed.  This  data  was  gathered  using  a  Particle  Image  Velocimetry  (PIV)  flow 
measurement  technique  with  a  trolling  motor  test  apparatus  in  the  test  section.  The  trolling 
motor  provides  similar  test  section  blockage  as  the  test  fixture  described  herein.  Note  that  the 
speed  reference  number  in  Figure  8  corresponds  to  the  output  frequency  from  the  impeller  motor 
drive  to  the  impeller  motor. 
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Flow  Speed  vs  Reference  Number 


Figure  8:  Test  Section  Flow  Speed  Determination 


Step  3  was  accomplished  by  operating  the  test  fixture  motor  drive  in  the  programmed  velocity 
mode  via  the  ASCII  command  line  of  the  Copley  Motion  Explorer  (CME)  software.  In  the 
programmed  velocity  mode,  a  rotation  speed  is  commanded  and  the  motor  drive  maintains  this 
speed  regardless  of  the  direction  of  energy  flow.  For  this  test  the  motor  is  acting  as  a  generator 
being  held  at  the  commanded  rotation  rate.  In  the  command  window  of  CME  it  was  observed 
that  the  RPM  was  being  held  to  the  commanded  RPM  +/-  2-3  RPM. 


Results  and  Comparison 

The  results  of  the  test  are  shown  in  Figure  9. 
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Figure  9  shows  the  following: 

1 .  There  is  good  agreement  between  predicted  and  experimental  data  for  tip  speed  ratios  (X) 
less  than  5. 

2.  On-design  predicted  perfonnance  almost  exactly  matches  the  experimental  on-design 
performance. 

3.  Experimental  and  predicted  performance  diverge  for  X  greater  than  5. 

As  a  result  of  the  experimental  results  shown  in  Figure  9,  OpenProp  is  being  revised  to  more 
accurately  predict  perfonnance  for  X  greater  than  ^Design-  It  is  thought  that  the  divergence  can  be 
accounted  for  by  implementing  a  more  sophisticated  model  of  bound  and  free  circulation  via 
lifting  surface  methods.  This  is  a  point  of  ongoing  work  in  OpenProp. 
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Chapter  4  -  Implementation  of  ABS  Steel  Vessel  Rules  for  Blade 
Thickness 


Manganese  Bronze 


Nickel- Manganese  Bronze 


Nickel- Aluminum  Bronze 


Manganese-Nickel- Aluminum-Bronze 

Stainless  Steel 


f,w-material  constants 


Upper  and  lower  blade 
surface  points  at  1 4  and  7/10 
radius 


P<> 70-pitch  at  0.7  radius 
divided  by  propeller  diameter, 
design  ahead  condition 


0*- 


Po  25-pitch  at  1 4  radius  divided 
by  propeller  diameter,  design 
ahead  condition 


a*-arca  of  expanded 
cylindrical  section  at 
0.25  Radius  mm*  or  in* 


V'OJ 


u , nr- 

modulus  coefficient 

it  r 

0.25 

Io-moment  of  inertia  of 
expanded  cylindrical  section 
at  '  4  radius  about  a  line 
passing  through  centroid 
parallel  to  pitch  line 

r  Up-Max.  nominal  distance  front  the 

o- 


moment  of  inertia  axis  to  points  of  the 
face  boundary  (tension  side)  in  mm  or  in 


C\  -  — Z-  section  area 
WT 

■•oefTicient  at  r-0.25 


c- =(i+i.5n„  x»r-s) 


_ 

W-Expanded  width  of 
cylindrical  section  at  1 4 
radius 

S' 

T-Max.  designed  thickness  of 
blade  section  at 1 4  radius 

H-Power  at  Rated  Speed 


5: 


f.w-material  constants 


N-Number  of  Blades 


R-RPM  at  rated  speed 


O 


_  S-Diameter  dependent 
factor 


D-Propeller  Diameter 


>«= 

4300*ra  Y  R  'i 

M  AlOoJ 

(3- 

1 

9 

a- expanded  blade  area 

divided  by  disc  area 

[Open Prop  does  not  allow  rake)  v. 

l0tSR  ~  S 

r  AH  ± 

im i 

''\C„C7W  * 

6- 


Figure  1:  Variable  Interrelationships  from  ABS  Steel  Vessel  Rules  for  Propellers 
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This  section  describes  the  implementation  of  the  American  Bureau  of  Shipping  (ABS)  steel 
vessel  rules  into  OpenProp  as  a  first  attempt  in  the  design  process  to  check  the  adequacy  of  the 
blade  dimensions  and  material  to  support  the  loads  they  will  carry.  The  output  of  OpenProp 
blade  structure  code  is  a  check  of  the  blade  thickness  at  the  quarter  span  section  against  the 
required  blade  thickness  at  the  quarter  span  section  as  determined  from  implementation  of  the 
steel  vessel  rules.  While  the  steel  vessel  rules  do  not  actually  calculate  a  stress  or  detennine  the 
operational  lifetime  of  a  propeller  they  do  take  these  quantities  into  consideration  as  evidenced 
by  the  rules  requirement  to  qualify  a  material  other  than  those  listed  for  service  in  a  classed 
vessel.  The  ABS  rules  also  represent  what  is  generally  required  in  order  to  class  a  vessel  with 
any  one  of  the  many  classification  societies  worldwide. 

Rule  Implementation  in  OpenProp 

The  OpenProp  module  which  implements  the  ABS  rules  for  propellers  does  so  in  a  way  which 
follows  the  flowchart  shown  in  the  figure  above.  User  input  for  this  module  is  only  the  material 
that  is  being  used  for  the  propeller  construction.  ABS  lists  five  different  materials  that  can  be 
used  for  propeller  manufacture;  these  are  listed  in  the  flowchart  above.  The  lines  of  code  which 
correspond  to  the  desired  material  must  be  uncommented  in  order  to  use  that  material  in  the 
calculations.  All  other  required  input  for  implementation  of  the  rules  is  automatically  extracted 
from  other  modules  of  OpenProp  or  calculated  within  the  blade  structure  module.  User  input  is 
highlighted  in  yellow;  input  from  other  modules  is  highlighted  in  green.  Since  other  OpenProp 
modules  use  the  SI  unit  system,  the  user  is  not  permitted  to  select  a  different  unit  system.  The 
output  of  the  structure  module  is  a  small  table  which  lists  the  section  thickness  at  the  quarter 
span  section  and  the  required  section  thickness  at  the  quarter  span  section,  as  calculated  from  the 
ABS  rules.  Propeller  redesign  is  necessary  if  the  required  blade  thickness  is  greater  than  the 
design  blade  thickness. 

Limitations 

In  its  current  version,  OpenProp  designs  fixed  pitch,  single  propellers  and  turbines  without  rake 
or  skew.  The  ABS  rules  for  propellers  allow  for  controllable  pitch,  rake  and  skew  but  the 
structure  module  developed  as  part  of  this  project  only  performs  the  calculations  for  fixed  pitch, 
single  propellers  without  rake  or  skew.  The  rules  used  to  develop  the  code  of  this  project  do  not 
cover  contra-rotating  propellers,  ducted  propellers  or  propellers  for  vessels  in  ice.  Additional 
structure  capability  could  easily  be  added  at  a  later  date  to  incorporate  the  ever  increasing 
capabilities  of  the  OpenProp. 

Moment  of  Inertia  Calculation 

The  bulk  of  the  code  to  implement  the  ABS  rules  for  propellers  is  used  to  detennine  the  moment 
of  inertia  of  the  designed  propeller  quarter  span  blade  section.  The  blade  structure  module  of 
OpenProp  imports  the  points  from  the  pressure  and  suction  sides  of  the  quarter  span  section.  All 
of  the  points  are  then  shifted  so  that  the  points  lie  in  the  first  quadrant  of  the  x-y  plane.  Shifting 
the  points  makes  the  determination  of  the  quarter  span  section  neutral  axis  easier.  The  code  then 
performs  a  trapezoidal  integration  for  the  pressure  and  suction  sides  separately  and  subtracts  the 
area  of  the  pressure  side  from  the  suction  side  so  that  only  the  enclosed  section  area  remains. 

The  moment  of  inertia  about  the  x-axis  is  then  calculated  and  the  parallel  axis  theorem  used  to 
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find  the  moment  of  inertia  about  the  neutral  axis.  A  flowchart  of  the  portion  of  the  code  which 
calculates  the  moment  of  inertia  is  shown  in  the  figure  below. 


Figure  2:  Code  Flowchart  to  Find  Section  Area  and  Moment  of  Inertia 
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Chapter  5  -  Calculation  of  Blade  Stress 


Theory 

A  relatively  simple  method  to  estimate  the  stress  on  a  propeller  blade  is  to  implement  beam 
bending  theory.  The  derivation  given  below  is  an  amplification  of  the  derivation  presented  in 
Kerwin  and  Hadler  (2010).  Kerwin  and  Hadler  also  include  some  historical  background  for  this 
method.  The  basic  assumptions  of  the  derivation  are: 

1 .  The  blade  acts  as  a  cantilevered  beam. 

2.  Axial  stresses  are  due  to  bending  and  centrifugal  forces. 

3.  Sheer  stresses  are  negligible. 

Figure  10  below  shows  a  propeller  blade  section  with  the  associated  inflow  velocities  and  lift 
force.  By  definition  the  lift  force,  dL,  is  always  perpendicular  to  the  total  inflow  velocity  V  .  dL 
is  responsible  for  both  thrust  and  torque  on  the  propeller  blades  and  propeller  shaft. 


Note  that  dL  is  always  perpendicular  to  V *  but  it  is  typically  not  perpendicular  to  the  chord  line. 
Therefore,  when  determining  the  component  of  dL  that  produces  thrust  and  the  component  of  the 
dL  which  produces  torque,  the  inflow  angle  ft  is  required,  not  the  blade  pitch  angle,  cpp.  The 
elemental  lift  at  a  blade  section  is  given  by  Equation  1 . 

dL  =  ^p(V*)2  CLcdr  (1) 


where 

dL  =  elemental  lift  on  a  blade  section 
p  =  fluid  density 

Cl  =  section  lift  coefficient  at  radius  r,  this  comes  from  the  lifting  line  calculation  in  OpenProp. 
c  =  section  chord  length  at  r 
dr  =  elemental  radial  span 
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Figure  11:  Blade  Section  Showing  Lift  Resolved  into  Axial  and  Tangential  Components 

From  Figure  1 1,  it  can  be  seen  that  the  axial  force,  Fa,  and  tangential  force,  Ft,  at  a  blade  section 
are  given  by: 


dFA  =  dL  cos(/?/  +  s )  (2) 
dFT  =  dL  sin(  fi,+e)  (3) 


where 

pi  =  inflow  angle 

e  =  Cd/Cl,  inflow  angle  correction  due  to  viscous  effects 
Cd  =  section  drag  coefficient 

Note  that  in  Figure  11,  the  point  of  application  of  dL  has  been  shifted  to  the  centroid  of  the 
section  and  is  no  longer  located  at  the  same  point  as  in  Figure  10.  This  is  done  to  simplify 
calculations.  dL  will  not  necessarily  be  located  at  the  section  centroid  but  at  a  point 
approximately  !4  of  the  distance  from  the  leading  edge  to  trailing  edge  on  the  chord  line  as 
shown  in  Figure  10.  The  fact  that  dL  does  not  act  through  the  section  centroid  means  that  dL  will 
produce  a  torque  about  the  span  line  of  the  blade.  This  torque  and  its  associated  sheer  stress  are 
assumed  to  be  negligible  along  with  all  other  shear  stresses. 

Both  Fa  and  Ft  produce  bending  moments  about  the  centroidal  axes.  Each  of  these  bending 
moments,  along  with  their  x  and  y  components,  is  shown  in  Figure  12.  The  equations  for  the 
bending  moment  are: 


dMA=(r-ro)dFA  (4) 
dMr  =  (r  —  ra)  dFT  (5) 


where 

rQ  =  section  radius  where  dM  is  being  calculated 
r  =  radius  of  section  producing  lift 
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Figure  12:  Bending  Moments  Components 

The  total  moments  produced  by  FA  and  Ft,  at  a  section  r0,  are  given  by: 

R  i 

ma  =  j~P(r*)2  CLc  cos(j8i+£)(r-ra)dr  (6) 

ro 

Mt  =  \^P(V*)2Clc  swift +e)(r-ra) dr  (7) 

ro 

Because  it  is  necessary  to  project  these  bending  moments  onto  the  centroidal  axes  of  the  section, 
blade  pitch  angle  is  required.  Projecting  the  total  bending  moments  onto  the  centroidal  axes,  the 
equations  become 


Mx  =Ma  cos((pP)+MT  sin(^p)  (8) 
MY  =Ma  sin((pp)-Mr  cos((pP)  (9) 
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Each  of  these  bending  moment  vectors  is  shown  in  Figure  13. 


Figure  13:  Total  Bending  Moments  about  Centroidal  Axes 

Additionally,  the  centrifugal  force  acting  at  each  section  contributes  to  the  overall  stress  at  the 
section.  The  elemental  centrifugal  force  acting  on  a  blade  from  an  adjoining  section  is  given  by 

dFc={2nnf1  r dm  (10) 

where 

dm  =  pbAdr  =  mass  of  blade  element 
Pb  =  propeller  blade  material  density 
A  =  section  area 
c  =  section  chord  length 
t  =  section  thickness 

Summing  the  contributions  of  all  adjoining  sections  to  the  Fc  at  the  section  of  interest,  the  total 
Fc  at  the  section  becomes 


R 

Fc  =  {2k  n)2  pbj  Ar  dr  (11) 

ro 

Since  the  blades  analyzed  using  the  above  method  do  not  contain  rake  or  skew,  which  would 
introduce  additional  bending  moments  from  Fc.  the  equation  for  the  stress  on  a  blade  section  can 
be  expressed  as: 


—  Mx  y  MYx  F 
a- - 2 - 2 — b  — 

lx.  Ij.  A 


(12) 
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Implementation 

The  following  paragraphs  are  an  explanation  of  the  method  used  to  calculate  blade  stresses.  The 
calculations  were  performed  on  a  propeller  which  was  designed  by  Epps  (2010a)  and  is  shown  in 
Figure  14. 


Figure  14:  Example  Propeller  Courtesy  of  Epps 

In  order  to  implement  the  equations  above  it  is  necessary  to  calculate  the  required  blade  section 
quantities.  Figure  15  and  Figure  17  illustrate  the  overall  procedure  for  determining  2D  blade 
section  area,  centroid  and  moments  of  inertia. 


Figure  15:  Distorted  Root  Section 


Figure  15  shows  the  visually  distorted  root  section  of  the  propeller  shown  in  Figure  14.  The 
section  is  distorted  for  illustrative  purposes;  the  undistorted  root  section  is  shown  in  Figure  16. 
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OpenProp  treats  the  blade  section  as  an  upper  and  lower  surface.  The  overall  procedure  for 
determining  blade  section  area  properties  consisted  of  determining  the  area  properties  of  the  area 
formed  by  the  upper  surface  and  the  x-axis  and  subtracting  the  area  properties  fonned  by  the 
lower  surface  and  the  x-axis.  This  subtraction  results  in  the  properties  of  the  enclosed  area 
shown  above. 


Figure  16:  Undistorted  Root  Section 


Figure  17:  Calculation  of  Elemental  Area  Properties 
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Figure  17  shows  a  characteristic  diagram  that  was  used  to  determine  elemental  area  properties 
which  were  summed  to  achieve  the  section  area  properties.  The  procedure  was: 

1 .  Determine  elemental  area 

2.  Calculate  elemental  centroid 

3.  Calculate  elemental  2nd  moment  of  area  about  both  x  and  y  axes 

4.  Sum  elemental  areas 

5.  Sum  2nd  moment  of  areas  about  x  and  y  axes 

6.  Calculate  section  centroid,  Equation  13 


14  ’  " 


7. 


Apply  parallel  axis  theorem  to  determine  2nd  moment  of  area  about  the  centroidal  axes, 
Equation  14. 


T  —T  —V  A 

1  Centroid  X  1  X  1  ^ total 


—  j  _  v1  A 

Centroid  Y  1  Y  ^ total 


(14) 


In  order  to  detennine  the  other  quantities  required  by  Equation  12,  the  integrals  were  turned  into 
discrete  sums  and  variables  from  the  propeller  design  were  used. 

Results 

The  results  of  the  analysis  performed  for  the  propeller  described  in  Epps  (2010a)  are  shown 
below  for  an  on-design  and  off-design  condition.  Figure  18  shows  the  stress  at  the  blade  root. 
As  expected,  the  blade  is  in  tension  on  the  pressure  side  and  compression  on  the  suction  side. 
Note  that  the  stresses  indicated  in  Figure  18  in  the  middle  of  the  root  section  are  interpolated 
stress.  Only  the  stresses  at  the  blade  surface  were  calculated  at  the  points  indicated. 
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In  Figure  19  through  Figure  22  tensile  stresses  are  considered  positive  and  compressive  stresses 
negative.  Figure  19  and  Figure  20  represent  the  on-design  condition  while  Figure  21  and  Figure 
22  represent  an  off-design  condition  as  specified  in  the  figure  titles.  As  expected,  the  off-design 
condition  chosen  shows  higher  stresses  than  the  on-design  condition  because  the  off-design 
condition  corresponds  to  a  point  where  the  propeller  is  producing  greater  thrust  and  torque. 
Greater  thrust  and  torque  results  in  higher  stress. 


Suction  Side  x10 

1 

0.5 

0 

■0.5 

-1 

-1.5 

Figure  19:  On-design  Suction  Side  Stress:  Js=0.75,  Vs=I.5m/s,  n=8rev/s,  D=0.25m 


Pressure  Side 


Figure  20:  On-design  Pressure  Side  Stress:  Js=0.75,  Vs=1.5m/s,  n=8rev/s,  D=0.25m 
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Figure  21:  Off-design  Suction  Side  Stress:  Js=0.40,  Vs=1.5m/s,  n=15rev/s,  D=0.25m 


Pressure  Side  x10 
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Figure  22:  Off-design  Pressure  Side  Stress:  Js=0.40,  Vs=1.5m/s,  n=15rev/s,  D=0.25m 


Carlton  (2007)  presents  isostress  contour  lines  taken  from  Finite  Element  Analysis  (FEA)  results 
for  various  propeller  types.  The  results  presented  above  agree  with  the  trends  presented  by 
Carlton  for  a  propeller  blade  without  skew.  Carlton  shows  highest  stress  near  the  blade  mid¬ 
chord  in  a  region  that  extends  close  to  the  tip  of  the  blade  and  a  decreasing  stress  as  one  moves 
away  from  the  mid-chord  to  the  blade  leading  and  trailing  edges. 
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Limitations 

There  are  several  assumptions  and  simplifications  that  are  made  in  the  method  discussed  above. 

Each  of  these  is  listed  below  with  suggestions  to  improve  the  calculation. 

1 .  Blade  Section  -  In  propeller  design  it  is  customary  to  present  the  blade  section  geometry 
as  the  unwrapped  section.  In  other  words,  the  blade  section  geometry  is  the  geometry 
one  would  obtain  if  the  curved  blade  section  were  laid  flat.  This  geometry  was  used  in 
the  calculation  of  the  blade  section  properties.  It  is  more  desirable  to  calculate  the  blade 
section  properties  for  blade  sections  that  were  taken  using  a  flat  cutting  plane  oriented 
perpendicular  to  the  span  line.  The  implementation  of  stress  calculations  would  be  more 
difficult  for  a  truly  flat  blade  section  since  OpenProp  does  not  develop  the  blade 
geometry  for  this  type  of  section  nor  is  the  hydrodynamic  data  valid  for  a  blade  section 
obtained  in  this  way.  The  OpenProp  geometry  code  could  be  rewritten  to  perform  an 
interpolation  in  order  to  obtain  the  points  to  create  a  truly  flat  blade  section  but  the 
problem  of  obtaining  blade  section  loading  remains.  Blade  loading  could  possibly  be 
obtained  using  the  cavitation  analysis  module  pressures  however  the  stress  code  would 
then  need  to  be  altered  to  be  something  similar  to  FEA. 

2.  Geometric  Property  Calculation  -  The  method  used  to  calculate  blade  section  properties 
was  essentially  a  trapezoid  rule  integration.  This  method  was  used  because  it  does  not 
limit  the  selection  of  the  number  of  points  used  to  design  the  blade,  however  other  more 
precise  methods  could  be  used. 

3.  Point  of  Load  Application  -  In  the  method  presented,  the  section  lift  force  is  applied  to 
the  section  centroid.  This  simplifies  the  code  because  the  actual  point  of  application  does 
not  have  to  be  determined  and  the  resulting  torsional  loads  can  be  ignored.  The  actual 
point  of  lift  force  application  could  be  approximated  as  the  14  chord  point  or  could  be 
calculted  more  precisely  by  analyzing  the  pressure  distribution  from  the  cavitation 
analysis  module.  If  the  load  point  is  corrected,  the  resulting  torsional  sheer  stresses  could 
be  calculated.  This  correction  is  probably  small  for  blades  without  skew. 

4.  Sheer  Stresses  -  All  sheer  stresses  on  the  blade  section  are  ignored.  The  stresses 
calculated  above  represent  bending  stress  from  lift  and  axial  stress  from  the  centrifugal 
force.  Sheer  stress  could  be  included  if  the  load  point  is  corrected  and  the  resulting 
torsional  stress  is  calculated  and  added  to  the  sheer  stress  resulting  from  the  pressure 
differential  between  the  blade  faces.  If  sheer  stresses  are  included,  the  principal  stresses 
should  be  calculated. 
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Chapter  6  -  Fatigue  Analysis 

By  definition  fatigue  failure  is  characterized  by  a  time  varying  load  whose  magnitude  is  smaller 
than  that  required  to  produce  failure  in  a  single  application,  Pook  (2007).  The  fatigue  analysis 
conducted  as  part  of  this  project  is  presented  in  two  steps. 

1 .  Identification  of  the  cyclic  loading 

2.  Application  of  a  fatigue  failure  theory. 

A  comprehensive  fatigue  analysis  is  characterized  by  many  subtleties  and  in  many  cases 
significant  experience  is  necessary  to  conduct  the  art  of  a  fatigue  analysis.  The  fatigue  analysis 
presented  here  is  intended  to  provide  a  method  by  which  a  fatigue  analysis  could  be  conducted 
on  a  propeller  or  turbine  at  the  beginning  of  the  design  process  to  ensure  the  estimated  fatigue 
life  meets  the  design  goal.  As  a  whole,  OpenProp  is  intended  to  be  a  design  tool  which  can  be 
used  to  provide  good  initial  propeller  and  turbine  designs.  As  additional  iterations  of  the  design 
process  are  completed  more  sophisticated  tools  for  propeller  design  will  become  necessary.  It  is 
in  this  spirit  of  providing  good  initial  design  estimates  that  the  fatigue  analysis  is  presented  here. 


Cyclic  Load 

For  a  propeller  or  turbine  the  source  of  the  varying  load  is  the  wake  that  it  operates  in.  Due  to 
the  presence  of  a  wake,  the  inflow  velocities  to  the  blades  are  not  uniform  in  magnitude  or 
direction.  As  a  blade  completes  a  revolution  it  will  pass  through  regions  of  various  velocity 
which  will  induce  varying  forces  on  the  blade.  A  propeller  will  typically  operate  in  a  wake  with 
greater  inflow  velocity  variation  than  a  turbine.  Because  a  propeller  operates  in  a  more  severe 
wake  environment  and  because  wake  data  is  more  readily  available  for  propellers  than  for 
turbines,  fatigue  analysis  for  a  propeller  was  performed. 

A  wake  for  a  single  screw  ship  is  shown  in  Figure  23.  This  wake  is  also  shown  in  Laskos  (2010) 
and  was  measured  by  Koronowicz,  Chaja  and  Szantyr  (2005).  This  figure  clearly  shows  a 
circumferential  variation  in  the  axial  inflow  velocity.  Typically  there  is  also  a  circumferential 
variation  in  the  tangential  inflow  velocity  but  this  variation  is  much  smaller  and  is  not  considered 
here.  This  is  shown  in  wake  profiles  of  Felli  and  Felice  (2005).  Figure  23  shows  the  ship  wake 
divided  into  twelve  sectors.  As  the  blade  passes  through  each  sector  it  is  assumed  to  fully 
develop  lift  commensurate  with  the  flow  velocity  in  that  sector.  This  assumption  makes  this 
analysis  a  quasi-steady  analysis.  In  each  sector,  the  circumferential  average  of  the  axial  inflow 
velocities  was  taken  at  the  same  radial  positions  that  were  used  in  the  propeller  design.  Each 
blade  section  is  subjected  to  a  different  inflow  velocity  which  results  in  a  different  C/_.  In  order 
to  determine  the  new  Cl  on  each  section  Equation  15  was  used. 

C  Lj- =C  Lo +27T  (Aoc)  (15) 


where 

Clo  =  original  lift  coefficient  in  the  design  condition 
Cl f  =  new  lift  coefficient  at  the  new  angle  of  attack 
Aa  =  change  in  angle  of  attack  from  design  condition 
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Figure  23:  Sectored,  Single  Screw  Ship  Wake 


Figure  24  and  Figure  25  show  how  a  change  in  the  axial  velocity  produces  a  change  in  the 
magnitude  and  direction  of  the  total  inflow  velocity.  The  analysis  also  assumes  that  u  a  and  u  t 
remain  constant. 
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Figure  26:  Pressure  Side  Blade  Stress  for  Each  Wake  Sector 
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Figure  26  above  shows  the  change  in  blade  stress  as  the  blade  passes  through  the  wake  sectors  of 
Figure  23.  As  expected,  the  highest  stresses  occur  in  sector  number  twelve  where  the  axial 
inflow  velocity  is  the  lowest.  The  lowest  axial  inflow  produces  the  largest  angle  of  attack  and 
lift  coefficient  and  subjects  the  blade  to  the  largest  amounts  of  lift  and  stress. 

Since  the  blade  stress  varies  considerably  across  the  blade  faces,  it  is  necessary  to  identify  the 
point  where  maximum  tensile  stress  occurs.  The  point  of  maximum  stress  for  this  propeller 
occurs  at  the  blade  root  at  the  point  identified  by  the  arrow  in  Figure  27. 


Figure  27:  Point  of  Maximum  Tensile  Stress 


In  Figure  28,  plots  of  the  maximum  blade  stress  versus  angular  blade  position  for  various  ship 
speeds  are  shown.  These  plots  also  identify  the  alternating  stress,  oa,  associated  with  each  blade 
stress.  Except  for  the  highest  ship  speeds,  <ra  is  relatively  low,  near  the  endurance  limit  for  nickel, 
aluminum  bronze,  as  shown  in  Figure  29. 
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Figure  28:  Maximum  Blade  Stress  versus  Angular  Position  for  Various  Ship  Speeds  Using  the  On-Design 

Advance  Coefficient  (Js=0. 75) 
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Fatigue  Failure 

Figure  29  shows  a  plot  of  <ra,  versus  number  of  reversals/cycles  to  failure  for  nickel,  aluminum 
bronze.  Data  for  this  figure  was  taken  from  Kerwin  and  Hadler  (2010);  detailed  alloy 
composition  and  test  condition  are  unknown.  Ideally,  one  would  design  a  propeller  such  that 
blade  stresses  were  minimized  in  order  to  increase  the  fatigue  life  of  the  propeller. 


When  performing  a  propeller  fatigue  analysis  it  is  critical  that  the  operational  profile  of  the  ship 
is  taken  into  consideration.  Figure  30  shows  an  operational  profile  for  a  warship  which  was 
taken  from  Gooding  (2009).  Since  the  propeller  analyzed  here  was  not  analyzed  for  such  a  wide 
spectrum  of  speeds,  Figure  31  was  used  in  the  example  calculation. 
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Example  Operational  Profile 


Figure  31:  Example  Operational  Profile  Used  for  Calculations 

With  the  assumptions  made  in  this  analysis,  there  is  a  direct  correlation  between  ship  speed  and 
blade  stress.  This  correlation  was  used  to  produce  Figure  32  below. 
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Figure  32:  Time  at  Various  Stress  Levels 


Miner’s  rule  was  used  to  predict  the  fatigue  life  of  the  propeller.  Miner’s  rule  is  simply  stated  in 
Equation  16. 


X-=i 

^  R; 


(16) 


where 

rt  =  actual  number  of  reversals  at  oa 

Rj  =  reversals  to  failure  at  oa ,  determined  from  Figure  29. 
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In  order  to  predict  the  fatigue  life,  additional  equations  are  necessary.  These  are  shown  below. 

ri=RPMitj  (17) 


where 

ti  =  time  spent  at  rotation  rate,  RPMj 

RPMj  =  rotation  rate  which  produces  desired  speed 


ti=xiT  (18) 


where 

Xj  =  fraction  of  total  time  spent  at  RPM; 
T  =  total  time  of  propeller  operation 


Substituting  Equation  17  and  Equation  18  into  Equation  16  and  solving  for  T,  one  obtains: 


T  = 


1 


RPM,  x , 


I 


R; 


(19) 


If  one  considers  the  blade  stress  at  speeds  below  25kts  to  be  of  infinite  life  then  the  fatigue  life  is 
180  days.  This  calculation  is  dominated  by  the  time  spent  at  30kts  which  is  probably  excessive 
when  comparing  Figure  3 1  and  Figure  30. 
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Limitations 

In  addition  to  the  limitations  discussed  at  the  end  of  Chapter  5,  there  are  some  additional 
limitations  that  are  specific  to  the  fatigue  analysis.  While  the  method  presented  here  is  useful  to 
obtain  an  initial  estimate  of  fatigue  life  there  are  many  more  factors  which  should  be  considered 
as  the  intial  estimate  is  refined. 

1 .  Fatigue  Data  -  Figure  29  shows  notional  fatigue  data  for  a  Ni-Al  Bronze.  It  is  unknown 
how  this  data  was  obtained  and  to  what  specific  alloy  it  applies.  It  would  be  desirable  to 
use  data  obtained  in  a  seawater  environment  for  the  specific  alloy  one  intended  to  use  to 
manufacture  a  propeller. 

2.  Miner’s  Rule  Coefficient  -  In  the  method  above,  it  was  assumed  that  when  the  sum  of 
Equation  16  reached  unity  that  the  material  would  fail.  This  is  a  reasonable  assumption 
but  other’s  have  found  that  the  coefficient  should  be  a  number  other  than  one  depending 
on  the  specific  material  type  (Sines  and  Waisman  1959).  Detennination  of  a  more 
accurate  value  for  this  coefficient  would  involve  an  extensive  test  program  that  should  be 
performed  for  more  refined  calculations. 

3.  Shaft  Reversals  -  Shaft  reversals  are  a  routine  operation  in  ship  maneuvering,  particularly 
when  the  ship  is  near  a  pier,  in  high  traffic  areas  or  navigating  waters  that  restrict  its 
turning  ability.  Kerwin  and  Hadler  (2010)  and  Carlton  (2007)  provide  a  limited 
discussion  on  the  effect  of  shaft  reversals  on  blade  fatigue  life.  A  shaft  reversal  can  be  a 
highly  stressing  event  for  the  propeller  blade  and  therefore  can  significantly  reduce  the 
blade  fatigue  life.  Calculating  the  blade  loads  during  a  shaft  reversal  would  require 
unsteady  analysis  and  is  beyond  the  scope  of  this  project  and  the  current  capabilities  of 
OpenProp. 

4.  Blade  Response  -  The  analysis  presented  here  assumes  that  the  blade  will  completely 
respond  to  the  flow  regime  of  each  sector.  It  is  unlikely  that  this  is  the  case.  The  analysis 
does  not  take  into  account  the  time  element  necessary  for  the  blade  to  fully  develop  its 
thrust  and  in  this  regard  is  a  conservative  estimate  of  fatigue  life. 


39 


Chapter  7  -  Test  Fixture  Design  and  Construction 

This  chapter  describes  the  design  and  construction  of  a  test  fixture  for  testing  propellers  and 
turbines.  The  test  fixture  described  in  this  chapter  was  specifically  designed  for  use  in  the 
hydrodynamics  laboratory  water  tunnel  at  MIT  but  can  also  be  used  in  a  tow  tank.  The 
limitations  of  the  test  fixture  are  given  in  the  Table  4. 


Limits 

Value 

Basis 

Torque 

6  ft-lbf 

Sensor  limitation 

Thrust 

50  lbf 

Sensor  limitation 

RPM 

1500  rpm 

Peak  capability  of  motor 

Current 

18  amps 

Peak  capability  of  motor 

Voltage 

240  V-AC 

300  V-DC 

Required  supply  voltage 

Maximum  controller  output  voltage 

Table  4:  Test  Fixture  Limitations 


The  design  philosophy  employed  for  this  test  fixture,  with  accompanying  justification  is  given 
below. 

1 .  Thrust  and  torque  sensor  must  be  the  limiting  component.  The  sensor  used  in  this  test 
fixture  is  on  loan  to  Professor  Richard  Kimball  from  the  US  Navy.  Searches  for  a 
commercially  available  sensor  capable  of  simultaneous  thrust  and  torque  measurement 
did  not  yield  any  devices  that  could  have  been  used  in  a  test  fixture  of  this  size.  Because 
of  the  limited  availability  of  useable  sensors,  it  was  decided  that  the  test  fixture  should  be 
limited  only  by  the  sensor. 

2.  Components  must  be  usable  in  other  test  fixtures.  Since  there  were  no  other  test  fixtures 
of  this  type  at  MIT,  this  design  constraint  meant  that  the  test  fixture  must  be  able  to  be 
disassembled  and  the  components  able  to  be  used  in  other  test  fixture  assemblies  that 
might  be  designed  by  students  in  the  future.  This  constraint  was  a  significant  driver  in 
the  selection  of  electrical  components,  manufacture  of  mechanical  components  and 
method  of  component  assembly. 

3.  Fixture  must  be  able  to  incorporate  a  high  resolution  encoder.  This  constraint  effected 
the  motor  and  encoder  selection  process. 

4.  Fixture  must  be  capable  of  use  in  both  a  tow  tank  and  water  tunnel.  This  constraint 
drives  the  maximum  allowable  overall  diameter,  length  and  standpipe  length  of  the  test 
fixture. 

Additional  details  concerning  how  the  design  philosophy  impacted  test  fixture  design  as  well  as 
the  final  test  fixture  configuration  are  given  in  the  sections  that  follow. 

MECHANICAL 

Thrust/Torque  Sensor 

The  thrust/torque  sensor  used  in  this  test  fixture  is  a  strain  gage  type  sensor.  The  sensor  uses  two 
sets  of  strain  gages;  one  set  to  measure  thrust  and  the  other  set  to  measure  torque.  The  strain 
gages  are  adhered  to  the  center  ring  shown  in  Figure  34,  which  is  covered  in  an  opaque  epoxy 
like  material.  The  presence  of  this  material  introduced  measurement  error  when  building  the 
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CAD  sensor  model  that  was  created  and  is  one  of  the  reasons  why  a  factor  of  safety  (FOS)  of  2 
was  used  when  determining  the  maximum  operational  torque  that  could  be  applied  to  the  sensor. 

Since  the  sensor  was  to  be  the  limiting  component,  it  was  necessary  to  characterize  the  thrust  and 
torque  capabilities  of  the  sensor.  In  order  to  determine  the  maximum  thrust  and  torque  that  the 
sensor  could  measure  without  damage,  a  determination  of  sensor  material  was  made  and  FEA  of 
the  sensor  was  performed.  For  the  purposes  of  this  test  fixture  design,  “damage”  is  defined  as  a 
load  condition  which  would  produce  yielding  in  the  sensor  material.  FEA  required  that  a  three 
dimensional  model  of  the  sensor  be  made,  this  model  is  shown  in  Figure  33. 


Figure  33:  Provided  Sensor 


While  measuring  the  sensor  to  determine  physical  dimensions  for  incorporation  into  the  model, 
an  inscription  of  501bf  was  found  on  one  end  of  the  sensor.  501bf  was  used  as  the  thrust  load  on 
the  sensor  in  the  FEA  analysis  in  order  to  determine  a  FOS.  The  result  of  the  FEA  showed  that 
the  sensor  can  withstand  a  501bf  axial  load  with  a  FOS  of  2.  The  calculated  stress  distribution 
resulting  from  a  501bf  load  is  shown  in  Figure  34. 


Figure  34:  Stress  from  Axial  Load  on  Sensor  (501bf  applied) 
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The  results  of  the  FEA  show  interference  between  the  center  ring  of  the  sensor  and  the  end  of  the 
sensor.  This  interference  is  a  result  of  the  large  scale  factor  necessary  to  make  the  sensor 
deflections  visible  and  does  not  represent  actual  interference  when  the  sensor  is  under  a  501bf 
thrust  load. 


In  order  to  determine  the  maximum  torque  that  the  sensor  could  carry,  a  separate  FEA  was 
conducted.  The  results  of  this  analysis  show  that  the  sensor  could  carry  12ft-lbf  without  damage. 
Application  of  a  FOS  of  2,  that  was  detennined  from  the  thrust  FEA,  limited  the  maximum 
torque  of  the  sensor  to  6ft-lbf.  A  FOS  of  2  is  reasonable  due  to  the  dimensional  error  present  in 
the  model  and  a  lack  of  validation  of  the  FEA  used  on  the  model  of  the  sensor.  A  picture  of  the 
stress  distribution  resulting  from  a  12ft-lbf  applied  torque  is  shown  in  Figure  35.  Note  that  the 
“handle”  that  is  present  in  the  picture  was  necessary  to  apply  a  torque  load  in  SolidWorks®  2007 
Education  Edition. 


Figure  35:  Stress  from  Torque  Load  on  Sensor  (12ft-lbf  applied) 


Output  Shaft  Configuration 

Three  options  were  considered  for  the  configuration  of  the  output  shaft  to  which  a  propeller  or 
turbine  could  be  attached  for  testing. 

1 .  A  tapered  shaft  capable  of  accepting  the  fittings  already  manufactured  and  located  in  the 
water  tunnel  laboratory. 

2.  A  straight  shaft  with  a  pin,  similar  to  that  used  for  propeller  attachment  to  trolling  motors. 

3.  A  straight  shaft  with  a  flat  side  machined. 

Option  1  was  undesirable  because  the  shaft  size  required  to  accommodate  the  taper  would  have 
required  larger  bearings  and  seals  for  the  shaft  which  would  have  increased  the  friction  resistance 
on  the  shaft  and  made  sealing  the  shaft  more  difficult.  Additionally,  a  larger  diameter  shaft  has 
greater  rotational  inertia  which  would  limit  the  rate  at  which  the  shaft  could  be  accelerated 
during  unsteady  tests. 
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Option  3  was  less  desirable  than  Option  2  because  of  the  complication  of  manufacturing 
propellers  with  a  set  screw  hole.  The  intended  manufacturing  technique  for  propellers  is  3D 
printing.  Propellers  manufactured  using  this  method  are  made  from  ABS  plastic.  Successfully 
creating  a  threaded  hole  into  this  material  with  sufficient  holding  power  for  a  set  screw  seemed 
unlikely.  A  second  problem  with  this  type  of  shaft  is  that  it  required  a  female  section  to  be  made 
in  the  propeller  hub  that  would  have  been  difficult  to  machine:  a  straight  cylindrical  hole  that 
changes  to  a  cylindrical  hole  with  a  flat.  Previous  experience  manufacturing  propellers  using  the 
3D  printing  technique  has  shown  that  it  is  difficult  to  achieve  a  hub  whose  outer  diameter  is 
concentric  to  the  drive  shaft  hole  outer  diameter.  Therefore  it  is  necessary  to  turn  the  propeller 
on  a  lathe  to  ensure  that  the  drive  shaft  will  easily  attach  to  the  propeller  with  minimal 
eccentricity  between  the  inner  and  outer  diameter  of  the  propeller  hub. 

Option  2  requires  that  every  propeller  have  a  slot  machined  in  the  hub  but  this  operation  is 
simple  using  an  end  mill  of  the  same  size  as  the  output  shaft  pin.  It  is  also  possible  to  print  the 
slot  in  the  hub  if  the  turbine  is  manufactured  using  a  rapid  prototyping  technique.  Option  2  also 
requires  that  the  end  of  the  drive  shaft  be  threaded  to  accept  a  nut  to  hold  the  propeller  against 
the  drive  shaft  pin,  however  these  are  external  threads  that  are  easy  to  manufacture.  For  these 
reasons,  Option  2  for  drive  shaft  configuration  was  chosen.  A  picture  of  the  shaft  is  given  Figure 
36. 


Drive  Shaft  Configuration 

The  test  fixture  design  described  in  this  paper  is  intended  to  be  used  to  test  both  propellers  and 
turbines.  Because  of  this  dual  use  capability,  it  is  necessary  that  the  fixture  be  able  to  measure 
and  support  axial  loads  in  two  directions.  Including  the  capability  to  support  axial  loads  in  two 
directions  also  protects  the  fixture  from  inadvertent  damage  should  a  load  be  applied  in  an  axial 
direction  for  which  the  fixture  was  not  designed. 
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Dual  axial  load  support  was  accomplished  by  using  two  tapered  roller  bearings  in  an 
arrangement  similar  to  the  front  wheel  bearing  assembly  on  an  older  automobile.  The  tapered 
roller  bearings  are  mounted  in  a  bearing  assembly  in  such  a  way  that  one  bearing  supports  the 
axial  load  in  one  direction  and  the  second  bearing  supports  the  axial  load  in  the  other  direction. 
The  drive  shaft  in  the  vicinity  of  these  bearings  is  threaded  and  slotted  to  accommodate  an  axle 
nut  and  star  washer.  The  nut  ensures  the  bearings  are  secured  in  the  bearing  housing  and  that  the 
axial  play  in  the  drive  shaft  can  be  adjusted.  The  slot  in  the  shaft,  in  combination  with  the  star 
washer,  ensures  that  the  nut  will  not  loosen.  A  picture  of  the  drive  shaft  and  bearing  assembly  is 
shown  in  Figure  37. 


Brush  Plate  O-Ring  Grooves 


Slip  Rings  and  Brushes  Sensor  Attachment  Flange 


Figure  37:  Driveshaft  and  Bearing  Assembly  with  Brush  Blocks  and  Slip  Rings 


The  smallest  diameter  on  the  driveshaft  was  determined  by  the  diameter  of  the  slip  ring  assembly. 
Due  to  the  shoulder  required  for  the  tapered  roller  bearings,  the  slip  ring  assembly  can  only  be 
installed  from  one  end  of  the  shaft.  The  end  of  the  shaft  over  which  the  slip  rings  must  be  moved 
to  reach  the  installation  location  was  made  slightly  smaller  than  the  slip  ring  diameter  in  order  to 
ease  slip  ring  installation.  The  drive  shaft  diameter  for  installation  of  the  slip  rings  is  only 
slightly  smaller  than  the  shaft  diameter  required  for  the  tapered  roller  bearings.  This  small 
change  in  diameter  meant  that  little  material  was  available  to  make  the  threads  for  the  axle  nut 
and  therefore  a  custom  nut,  washer  and  thread  configuration  had  to  be  manufactured. 


Housings 

Two  assembly  methods  were  considered  for  the  external  housings. 

1 .  Threaded  assembly 

2.  Shoulder  fasteners 
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Using  a  threaded  assembly  has  the  advantage  of  minimizing  the  number  of  water  leakage  paths 
into  the  fixture  and  the  number  of  o-rings  required  during  assembly.  The  problem  with  a 
threaded  assembly  is  that  the  threads  can  be  difficult  to  manufacture,  particularly  for  internal 
threads  that  run  deep  into  the  part,  and  large  diameter  threads  are  prone  to  seize  in  stainless  steel. 
The  problem  with  the  shoulder  fastener  assembly  method  is  that  the  number  of  leakage  paths  and 
o-rings  required  is  significant  and  assembly  requires  that  the  components  be  precisely  positioned 
prior  to  the  installation  of  the  shoulder  fasteners.  The  shoulder  fastener  assembly  method  was 
chosen  for  ease  of  manufacture  and  the  problem  of  water  leakage  paths  was  mitigated  by 
installing  the  shoulder  fasteners  between  a  set  of  o-rings  on  the  housing  diameters. 


ELECTRICAL 


Slip  Rings 

The  torque/thrust  sensor  uses  two  sets  of  strain  gages  for  load  measurement.  There  are  several 
possible  methods  to  transmit  the  data  signal  from  the  sensor  for  recording.  The  method  chosen 
for  this  fixture  was  to  amplify  the  data  signal  at  the  sensor  and  then  use  a  set  of  slip  rings  and 
brushes  to  conduct  this  signal  to  a  point  where  a  data  acquisition  system  could  be  attached.  This 
method  was  chosen  for  its  simplicity. 


The  details  of  the  slip  ring  assembly  are  included  in  the  appendix.  Six  slip  rings  are  required  for 
sensor  operation.  Two  are  necessary  to  power  both  sets  of  strain  gages,  four  rings  are  necessary 
for  data  signal  transmission.  The  slip  assembly  used  in  this  test  fixture  has  eight  slips  rings  in 
order  to  allow  for  future  growth  and  to  provide  alternate  slip  rings  should  some  become  unusable. 
Each  slip  ring  has  four  brushes  riding  on  it,  two  from  each  brush  block.  The  brushes  from  each 
brush  block  are  soldered  together  so  that  four  brushes  are  connected  to  each  slip  ring.  Four 
brushes  per  slip  ring  are  used  in  order  to  minimize  the  electrical  resistance  between  the  brushes 
and  the  slip  rings.  A  photograph  of  the  slip  rings  and  brushes  installed  in  the  test  fixture  is 
shown  in  Figure  38. 


Brush 
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High  Vacuum  Grease 


Electrical  Tape  and  Zip  Tie 


Leakage  Alarm  Sensor 


Figure  38:  Installed  Slip  Rings  and  Brushes 
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In  addition  to  the  slip  rings  and  brushes,  Figure  38  shows  two  red  wires  near  the  bottom  of  the 
photograph.  These  wires  are  part  of  a  leakage  alann  system.  If  the  ends  of  the  two  red  wires  are 
shorted,  an  audible  alarm  sounds  indicating  leakage  into  the  test  fixture.  The  yellow  electrical 
tape  and  zip  tie  are  present  to  secure  the  slip  ring  and  strain  gage  amplifier  wiring  to  the 
driveshaft.  The  connection  of  the  slip  ring  wires  and  the  stain  gage  amplifiers  was  made  outside 
of  the  hollow  drive  shaft  in  order  to  ease  assembly.  The  white  substance  on  the  end  of  the  tube  is 
a  high  vacuum  silicone  based  grease  that  is  applied  to  the  surfaces  prior  to  assembly  in  order 
ease  assembly  and  as  an  additional  measure  to  prevent  leakage. 


Amplifiers 

Inside  the  thrust/torque  sensor  are  two  amplifiers.  One  amplifier  is  for  the  thrust  data  signal  and 
the  other  is  for  the  torque  data  signal.  These  amplifiers  are  mounted  inside  a  piece  of  foam 
which  is  pressed  into  the  sensor.  The  amplifiers  that  were  purchased  are  designed  for  strain  gage 
signal  amplification  for  the  motor  sports  industry  and  therefore  represent  a  rugged  option  for 
signal  amplification.  Data  signal  amplification  takes  place  as  close  to  the  sensor  as  possible  in 
order  to  limit  the  data  signal  transmission  loss  and  to  prevent  the  signal  to  noise  ratio  of  the  data 
signal  from  becoming  too  low  for  practical  use.  Additional  amplifier  details  are  included  in  the 
appendix. 


Motor 

The  desire  to  use  the  test  fixture  in  the  water  tunnel  limits  the  maximum  allowable  diameter  of 
the  test  fixture.  Previous  experience  with  trolling  motors  in  the  water  tunnel  yielded  good  results. 
Trolling  motor  diameters  ranged  from  3.5  to  4  inches;  therefore  the  maximum  allowable  test 
fixture  diameter  was  set  to  4  inches.  A  maximum  diameter  of  4  inches  significantly  restricts  the 
available  options  for  motor  selection. 

Another  consideration  in  motor  selection  is  the  ability  of  the  motor  to  also  act  as  a  generator  in 
order  to  serve  as  a  load  for  turbine  testing.  The  requirement  to  also  act  as  a  generator  further 
limits  the  choice  of  motor  to  those  of  a  pennanent  magnet  design.  Although  it  is  possible  to  use 
a  motor  without  permanent  magnets  installed,  the  complication  arising  from  supplying  both  the 
stator  and  rotor  with  electric  current  was  deemed  excessive  for  a  test  fixture. 

The  motor  selected  for  this  test  fixture  is  a  Parker  kit  motor,  K089300.  This  motor  is  a  DC 
brushless  motor;  the  specific  model  selected  also  contains  integral  commutation.  In  selecting  a 
motor,  it  was  desirable  to  select  a  motor  such  that  the  sensor  remained  the  limiting  component  in 
the  design.  Therefore  a  motor  capable  of  torque  in  excess  of  16N-m  (12ft-lbf)  was  selected.  A 
test  fixture  using  a  standard  motor  required  a  test  fixture  6  inches  in  diameter. 

A  standard  motor  with  the  desired  torque  speed  characteristics  exceeded  the  maximum  allowable 
diameter  because  a  standard  motor  comes  with  a  face  plate  on  one  end  and  electrical  connectors 
on  the  other.  The  standard  motor  would  have  required  customization  to  remove  the  electrical 
connectors  and  change  the  mounting  configuration  to  a  face  frame  mount.  In  addition  to  the 
complication  and  cost  of  performing  the  customization,  supplying  the  motor  with  current  would 
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have  also  been  challenging  because  the  wires  would  have  had  to  pass  by/through  the  face  plate 
area  in  order  to  be  routed  to  the  standpipe  for  passage  out  of  the  test  fixture  to  the  electrical 
supply. 

The  K089  series  motor  has  a  maximum  diameter  of  3.5  inches  which  allows  the  motor  to  be 
attached  inside  a  tube  with  a  maximum  outside  diameter  of  4  inches.  The  K089300  was  selected 
because  it  is  the  highest  torque  motor  listed  in  the  catalogue  for  this  series.  Parker  frameless  kit 
motors  generate  additional  torque  in  a  given  series  by  increasing  the  length  of  the  stator  windings 
and  rotor.  A  kit  motor  has  the  following  advantages: 

1 .  Able  to  fit  in  a  tighter  package 

2.  Use  of  one  shaft  for  both  the  motor  and  driveshaft  which  eliminates  the  need  for  a  shaft 
coupling  between  motor  shaft  and  propeller  drive  shaft 

3.  Stator  windings  are  in  direct  contact  with  the  test  fixture  housing  which  allows  efficient  heat 
transfer  out  of  the  stator  windings  and  into  the  fluid  surrounding  the  test  fixture. 

The  disadvantages  of  using  a  kit  motor  are: 

1 .  Attachment  of  the  motor  into  the  fixture  required  that  additional  holes  had  to  be  drilled 
into  the  motor  housing  which  meant  an  increase  in  the  probability  of  a  leak  into  the  fixture. 

2.  Kit  motors  do  not  come  with  a  high  resolution  angle  encoder  installed  like  the  standard  off 
the  shelf  motors. 

The  calculated  torque  speed  curve  for  the  K089300  is  given  in  Figure  39.  In  both  figures  the 
dashed  lines  represent  continuous  operation  while  the  solid  lines  represent  peak  or  intermittent 
operation.  The  linear  negative  slope  in  the  torque-speed  curve  is  based  on  preventing  the  motor 
windings  from  overheating  due  to  excess  current. 
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Figure  39:  K089300-7Y  Torque  Speed  Curve 
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0. 


Output  Power 


Figure  40:  Output  Power  Capability 

The  current  limitation  of  the  motor  is  shown  in  Figure  41. 


Controller 

In  selecting  a  controller  it  was  desirable  to  select  a  controller  which  could  serve  as  both  a  motor 
controller  and  load  controller,  allowed  for  test  fixture  growth  and  had  limited  EMI  to  prevent 
noise  in  the  data  signal.  Students  at  the  University  of  Maine  have  built  and  used  a  test  fixture  for 
cross  flow  turbines  that  used  a  Copley  Xenus  XTL-230-40  controller.  They  have  been  very 
pleased  with  the  overall  performance  of  their  test  fixture,  particularly  the  low  electrical  noise 
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generated  by  the  controller.  For  these  reasons  the  same  controller  and  electrical  layout  were 
selected  for  this  test  fixture.  A  photograph  of  the  electrical  components  in  the  enclosure  is 
shown  in  Figure  42. 
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Figure  42:  Electrical  Components 
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Figure  43:  Schematic  of  Enclosure  Electrical  Components 


A  schematic  of  the  components  that  are  located  inside  the  enclosure  is  shown  in  Figure  43.  Note 
that  all  connection  points  to  the  motor  drive  are  not  shown,  only  those  connections  that  are  used 
are  shown.  Strain  gage  wiring  is  a  different  circuit  and  is  not  shown  in  Figure  43. 

CONSTRUCTION 

The  material  used  in  the  construction  of  the  parts  of  the  test  fixture  that  could  be  in  contact  with 
the  water  is  stainless  steel.  Depending  on  the  part,  the  alloy  is  either  a  303  or  304  stainless  steel. 
These  alloys  were  selected  for  their  combination  of  corrosion  resistance  and  machinability. 

Their  corrosion  resistance  will  be  sufficient  for  use  in  a  fresh  water  environment,  however 
prolonged  use  in  chlorinated  water  and  use  in  saltwater  should  be  avoided  to  prevent  corrosion. 
All  components  in  the  test  fixture  are  non-magnetic  with  the  exception  of  the  drive  shaft  and 
propeller  shaft  which  became  slightly  magnetic  as  a  result  of  the  machining  process.  A 
photograph  of  the  completed  test  fixture  in  operation  during  turbine  testing  are  shown  in  Figure 
44. 
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Figure  44:  Completed  Test  Fixture  in  Operation 
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Chapter  8  -  Conclusions  and  Recommendations 
Conclusions 

In  this  thesis  a  propeller  and  turbine  test  fixture  was  designed,  constructed  and  used  to  test  a 
model  scale  axial  flow  hydrokinetic  turbine.  Additionally,  computer  codes  were  written  to 
implement  ABS  Steel  Vessel  Rules  for  propellers,  calculate  blade  stresses  and  perform  an  initial 
estimate  of  fatigue  life  for  a  propeller  blade. 

Results  from  turbine  tests  show  that  the  test  fixture  produces  data  signals  with  very  little 
interference.  Test  results  also  show  that  the  OpenProp  perfonnance  predictions  for  the  turbine 
design  point  are  very  good  but  that  additional  work  is  necessary  to  improve  the  off-design 
performance  predictions  of  OpenProp. 

Plots  of  blade  stress  distributions  obtained  using  the  methods  described  above  show  good 
agreement  with  those  of  Carlton  (2007)  which  were  obtained  using  more  sophisticated 
techniques.  The  method  used  to  calculate  blade  stress  was  easily  extended  to  perform  a  fatigue 
analysis  to  obtain  an  initial  estimate  of  a  propeller  blade  turning  in  a  wake. 


Recommendations  for  Further  Work 

As  currently  constructed,  the  electrical  enclosure  for  the  test  fixture  must  be  operated  with  the 
door  open  in  order  to  pass  the  power  and  control  cables.  The  cables  should  be  connectorized  and 
outlets  attached  to  the  enclosure  so  that  the  door  can  be  shut  during  testing. 

The  test  procedure  for  collecting  performance  data  is  a  manual  process  using  LabView®. 

Lab  View®  has  the  capability  to  automate  the  process  in  a  way  that  would  make  data  collection 
much  quicker  and  the  results  more  reliable.  The  test  fixture  could  be  controlled  using  the  CME2, 
ASCII  command  interface  that  is  part  of  the  motor  controller. 

Unsteady  tests  would  require  a  high  resolution  shaft  position  encoder.  The  test  fixture  was 
designed  to  allow  for  the  inclusion  of  a  Renishaw™  encoder.  Two  additional  pieces  would  have 
to  be  manufactured;  drawings  for  these  parts  are  included  in  the  appendix. 

In  order  to  implement  a  method  to  calculate  blade  stress,  several  simplifications  were  made. 

More  precise  results  could  be  obtained  by  implementing  an  FEA  approach.  OpenProp  already 
calculates  all  of  the  necessary  geometry  and  load  data  that  could  be  used  in  an  FEA  code.  This 
method  should  be  implemented  with  default  settings  that  are  tailored  to  the  geometry  of  propeller 
and  turbine  blades  such  that  a  person  with  little  FEA  experience  could  obtain  valid  results. 

The  fatigue  results  presented  assumed  a  quasi-steady  blade  response  as  the  blade  passed  from 
one  wake  sector  to  another.  More  sophisticated  techniques  exist  to  estimate  blade  load  as  it 
passes  into  regions  of  various  flow  speed  in  the  wake.  These  techniques  should  be  implemented 
in  order  to  refine  the  fatigue  analysis  predictions. 
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Appendix  A  -  Codes 

Moment  of  Inertia  Calculation 

function  [Mpl,  Ixc,  Iyc,  Ixyc,  A,  Xbar,  Ybar,  xl,  yl,  xu,  yu]  =  MomentofInertia(xl,xu,yl,yu) 

[Mp  1  ,Np]  =  size(xu); 

%  Calculation  of  Section  Area  and  Centroid 
for  m=l:Mpl 

yshift  =  abs(min(yl(m,:)));  %Distance  to  shift  all  y  points  so  that  all  are  positive 
yu(m,:)  =  yu(m,:)  +  yshift;  %Shift  of  upper  surface  y  points 

yl(m,:)  =  yl(m,:)  +  yshift;  %Shift  of  lower  surface  y  points 

xshift  =  abs(min(min(xu(m,:)),min(xl(m, :))));  %Distance  to  shift  all  x  points  so  that  all  are 

positive 

xu(m,:)  =  xu(m,:)  +  xshift;  %Shift  of  upper  surface  x  points 

xl(m,:)  =  xl(m,:)  +  xshift;  %Shift  of  lower  surface  x  points 

end 

dxu  =  abs(diff(xu,l,2)); 
dxl  =  abs(diff(xl,  1 ,2)); 
dyu  =  diff(yu,l,2); 
dyl  =  diff(yl,l,2); 

Ybar  =  zeros(  1  ,Mp  1 ); 

Xbar  =  Ybar; 

Ixc  =  Ybar; 

Iyc  =  Ybar; 

A  =  Ybar; 

Ixyc  =  Ybar; 

for  m=l:Mpl 

hru  =  zeros(  1  ,(Np- 1 )); 
hrl  =  hru; 

htu  =  zeros(l,(Np-l)); 
htl  =  htu; 

xctu  =  zeros(l,(Np-l)); 
xctl  =  xctu; 

for  n=l:(Np-l) 

hru(n)=min(yu(m,n),yu(m,n+l));  %Height  of  upper  surface  elemental  rectangle 
htu(n)=max(yu(m,n),yu(m,n+ 1));  %Height  of  upper  surface  elemental  trapezoid 
hrl(n)=min(yl(m,n),yl(m,n+l));  %Height  of  lower  surface  elemental  rectangle 
htl(n)=max(yl(m,n),yl(m,n+l));  %Height  of  lower  surface  elemental  trapezoid 
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if  dyu(m,n)<0 

xctu(n)  =  xu(m,n)  +  2*dxu(m,n)/3; 
else 

xctu(n)  =  xu(m,n)  +  dxu(m,n)/3; 
end 

if  dyl(m,n)>0 

xctl(n)  =  xl(m,n)  +  2*dxl(m,n)/3; 
else 

xctl(n)  =  xl(m,n)  +  dxl(m,n)/3; 
end 
end 

xcru  =  xu(m,l:(Np-l))+dxu(m,:)/2; 
xcrl  =  xl(m,l:(Np-l))+dxl(m,:)/2; 


%Distance  from  y-axis  to  upper  surface 
elemental  triangle  centroid 

%Note:  Value  depends  on  whether  left  or  right 
side  of  triangle  is  higher 


%Distance  from  y-axis  to  lower  surface 
elemental  triangle  centroid 

%Note:  Value  depends  on  whether  left  or  right 
side  of  triangle  is  higher 


%Distance  from  y-axis  to  upper  surface 
elemental  rectangle 

%Distance  from  y-axis  to  lower  surface 
elemental  rectangle 


aru  =  dxu(m,:).*hru; 

atu  =  dxu(m,:).*(htu-hru)/2; 

arl  =  dxl(m,:).*hrl; 

atl  =  dxl(m,:).*(htl-hrl)/2; 


%Elemental  upper  surface  rectangle  area 
%Elemental  upper  surface  triangle  area 
%Elemental  lower  surface  rectangle  area 
%Elemental  lower  surface  triangle  area 


ycru  =  hru/2; 

ycrl  =  hrl/2; 

yctu  =  hru+(htu-hru)/3; 

yeti  =  hrl+(htl-hrl)/3; 


%Distance  from  x-axis  to  upper  surface  elemental  rectangle 
centroid 

%Distance  from  x-axis  to  lower  surface  elemental  rectangle 
centroid 

%Distance  from  x-axis  to  upper  surface  elemental  triangle 
centroid 

%Distance  from  x-axis  to  lower  surface  elemental  triangle 
centroid 


Mxsu  =  sum(ycru.*aru  +  yctu.*atu);  %lst  moment  of  upper  surface  about  x  axis 
Mxsl  =  sum(ycrl.*arl  +  yctl.*atl);  %lst  moment  of  lower  surface  about  x  axis 

Mxs  =  Mxsu  -  Mxsl; 


Mysu  =  sum(xcru.*aru  +  xctu.*atu);  %  1  st  moment  of  upper  surface  about  y  axis 
Mysl  =  sum(xcrl.*arl  +  xctl.*atl);  %lst  moment  of  lower  surface  about  y  axis 
Mys  =  Mysu  -  Mysl; 
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Au  =  sum(aru  +  atu);  %Area  of  upper  surface  (x  axis  to  upper  surface) 

A1  =  sum(arl  +  atl);  %Area  of  lower  surface  (x  axis  to  lower  surface) 

A(m)  =  Au  -  Al; 

Ybar(m)  =  Mxs/A(m);  %Distance  to  centroid  from  x-axis 

Xbar(m)  =  Mys/A(m);  %Distance  to  centroid  from  y-axis 

%Uncomment  lines  below  to  see  a  section  graph  with  centroidal  axes 
%figure(m) 

%plot(xu(m,:),yu(m,:),xl(m,:),yl(m,:),'b') 

%line([min(xu(m,:)),max(xu(m,:))],[Ybar(m),Ybar(m)],'ColorVr','LineWidth',2,'LineStyleV- 

-’) 

%line([Xbar(m),Xbar(m)],[min(yl(m,:)),max(yu(m,:))],'ColorVr','LineWidth',2,'LineStyleV- 

-’) 

%axis  equal 
%grid  on 

%  Calculation  of  Section  Moment  of  Inertia 

ixru  =  dxu(m,:).*hru.A3/3; 

ixyru  =  aru.*ycru.*xcru; 

ixtu  =  dxu(m,:).*(htu-hru).A3/36  +  atu.*yctu.A2; 

ixytu  =  atu.*yctu.*xctu; 

ixrl  =  dxl(m,:).*hrl.A3/3; 
ixyrl  =  arl.*ycrl.*xcrl; 

ixtl  =  dxl(m,:).*(htl-hrl).A3/36  +  atl.*yctl.A2; 
ixytl  =  atl.*yctl.*xctl; 

iyru  =  hru.*dxu(m,:).A3/12  +  aru.*xcru.A2; 
iytu  =  (htu  -  hru).*dxu(m,:).A3/36  +  atu.*xctu.A2; 
iyrl  =  hrl.*dxl(m,:).A3/12  +  arl.*xcrl.A2; 
iytl  =  (htl  -  hrl).*dxl(m,:).A3/36  +  atl.*xctl.A2; 

lx  =  sum(ixru  +  ixtu)  -  sum(ixrl  +  ixtl); 

Iy  =  sum(iyru  +  iytu)  -  sum(iyrl  +  iytl); 

Ixy  =  sum(ixyru  +  ixytu)  -  sum(ixyrl  +  ixytl); 

Ixc(m)  =  lx  -  A(m)*Ybar(m)A2; 

Iyc(m)  =  Iy  -  A(m)*Xbar(m)A2; 

Ixyc(m)  =  Ixy  -  A(m)*Xbar(m)*Ybar(m); 
end 
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Centrifugal  Force  Calculation 

function  [omega,  Fc,  Rdif]  =  CentrifugalForce(Mpl,  N,  R,  A,  RC,  DR) 

%Initialization  of  Variables 
Fc  =  zeros(l,Mpl-l); 

Rdif  =  zeros(Mp  1  - 1  ,Mp  1  - 1 ); 

R1  =  Rdif; 

omega  =  2*pi*N/60;  %[rev/s]  Rotation  rate 

gamma  =  1024. 16;  %[kg/mA3]  Density  of  ABS  plastic 

%  8200;  %[kg/mA3]  Approximate  density  of  Ni-Al-Bronze 

M  =  (omega)A2  *  gamma  *  RA2;  %Multiplier  used  below 

%Calculation  of  centrifugal  force  on  each  section 
for  m  =  l:(Mpl-l) 

Fc(m)  =  M  *  sum(A(m:end-l).*RC(m:end).*DR(m:end));  %[N] 
Rl(m,:)  =  RC  -  RC(m); 

I  =  find(Rl(m,:)>0); 

Rdif(m,I)  =  Rl(m,I); 
end 

end 


58 


Stress  Calculation 


function  [s]  =  Stress(CD,  CL,  BetalC,  Mpl,  Np,  xu,  yu,  xl,  yl,  rho,  Vs,  R,  VSTAR,  CoD, 
Rdif,  DR,  theta,  Xbar,  Ybar,  Ixc,  Iyc,  Ixyc,  Fc,  A) 

%Uncomment  lines  below  to  use  Kerwin  and  Hadler  method  exactly 
%  eps  =  CD./CL; 

%  S  =  sin(BetalC-eps); 

%  C  =  cos(BetalC-eps); 

S  =  sin(BetalC);  %Factors  for  use  below 
C  =  cos(BetalC);  %Factors  for  use  below 

%Initialization  of  Variables 
INTQ  =  zeros(Mpl-l,Mpl-l); 

INTT  =  INTQ; 

MQ  =  zeros(l,Mpl-l); 

MT  =  MQ; 

Mxo  =  MQ; 

Myo  =  MQ; 

s  =  zeros(Mpl-l,2*Np); 

%Concatenation  of  upper  and  lower  section  curves  into  a  single  curve 

Xu  =  xu 
Yu  =  yu(:,:); 
xs  =  cat(2,Xu,fliplr(xl)); 
ys  =  cat(2,Yu,fliplr(yl)); 

%Uncomment  lines  below  to  see  a  plot  of  root  blade  section 
%  axes('fontweighf  ,’bold') 

%  hold  on 

%  plot(xs(  1  ,:),ys(  1  ,:),'-ks','LineWidth',2) 

%  axis  equal 

%  title(’Root  Section  Plof,’Fontweighf , 'bold',  ’Fontsize’,  14) 

%  xlabel('X  (m)’,’Fontweighf  ,’bold’, ’Fontsize’,  12) 

%  ylabel(’Y  (m)’,’Fontweighf  ,’bold', ’Fontsize’,  12) 

M  =  rho*VsA2*RA3;  %Multiplier  used  below 

for  m=l:(Mpl-l) 

%Uncomment  two  lines  below  to  use  Kerwin  and  Hadler  method  exactly 

%INTQ(m,:)  =  M*  (VSTAR  A2  .*  CL  *  CoD  .*  Rdif(m,:)  AS  ADR); 

%INTT(m,:)  =  M*  (VSTAR. A2  *  CL  A  CoD  .*  Rdif(m,:)  AC  ADR); 

INTQ(m,:)  =  M*  (VSTAR. A2  A  (CL  AS  +  CDAC)A  CoD  A  Rdif(m,:)ADR); 
INTT(m,:)  =  M*  (VSTARA2  A  (CL  AC  -  CDAS)A  CoD  A  Rdif(m,:)ADR); 

MQ(m)  =  sum(INTQ(m,:)); 
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MT(m)  =  sum(INTT(m,:)); 

Mxo(m)  =  MT(m)*cos(theta(m))  +  MQ(m)*sin(theta(m)); 

Myo(m)  =  MT(m)*sin(theta(m))  -  MQ(m)*cos(theta(m)); 

xsdiff(m,:)  =  xs(m,:)  -  Xbar(m); 

ysdiff(m,:)  =  ys(m,:)  -  Ybar(m); 

s(m,:)  =  ((-Mxo(m)*Iyc(m)  +  Myo(m)*Ixyc(m))*ysdiff(m,:)  -  (-Mxo(m)*Ixyc(m)  + 
Myo(m)*Ixc(m))*xsdiff(m,:))  /  (Ixc(m)*Iyc(m)  -  Ixyc(m)A2)  +  Fc(m)/A(m); 

%Uncomment  line  below  to  use  a  more  exact  equation  for  stress  which 
%takes  into  account  the  product  of  inertia 

%s(m,:)  =  ((-Mxo(m)*Iyc(m)  +  Myo(m)*Ixyc(m))*ys(m,:)  -  (-Mxo(m)*Ixyc(m)  + 
Myo(m)*Ixc(m))*xs(m,:))  /  (Ixc(m)*Iyc(m)  -  Ixyc(m)A2)  +  Fc(m)/A(m); 

%  Uncomment  lines  below  for  plots  of  stress  on  each  blade  section 
%  figure(m) 

%  %  plot3(xs(m,:),ys(m,:),s(m,:),'rs') 

%  %  grid  on 

%  %  xlim([min(xs(m,:)),max(xs(m,:))]) 

%  %  ylim(xlim) 

%  patch(xs(m,:),ys(m,:),s(m,:)) 

%  colonnap(jet) 

%  colorbar 
%  grid  on 
%  axis  equal 
end 

clear  CL 
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Blade  Stress  Plots 


function  [  ]  =  Plot_Blade_Contours(X3D,Y3D,Z3D,s,plottitle) 

[rows,cols]=size(X3D);  %X3D  is  from  the  geometry .m  module 
Mp  =  rows- 1 ;  %Number  of  blade  sections 

%oConcatenate  matrices  to  create  vertex  matrix  for  patch  function.  This  is  necessary  because 
the  patch  function  expects  a  matrix  whose  rows  are  the  vertices  in  x,y,z  coordinates. 

colx  =  X3D(1,:)'; 
coly  =  Y3D(1,:)'; 
colz  =  Z3D(1,:)'; 
colS  =  s(  1 
for  n=2:rows-l 

colx  =  vertcat(colx,X3D(n,:)'); 
coly  =  vertcat(coly,Y3D(n,:)'); 
colz  =  vertcat(colz,Z3D(n,:)'); 
colS  =  vertcat(colS,s(n,:)'); 

end 

%Create  face  matrix  for  patch  function-  this  tells  the  patch  function  how  to  connect  the 
vertices  to  create  a  face.  This  code  uses  a  square/rectangular  face. 

F(:,l)  =  l:cols*(Mp-l); 

F(:,2)  =  F(:,l)  +  1; 

F(:,3)  =  F(:,l)  +  (cols+1); 

F(:,4)  =  F(:,l)  +  cols; 

%Create  special  rows  -  the  pattern  of  face  vertices  "wraps, "  these  lines  make  the  pattern  wrap 
properly 

m  =  0:cols:cols*(Mp-l); 
for  n  =  1  :length(m) 
if  m(n)>l 

F(m(n),l:4)  =  [F(m(n),l),  F(m(n),l)-(cols-l),  F(m(n),l)+1,  F(m(n),l)+4]; 
end 

end 

%>Remove  extra  face  matrix  rows  -  This  removes  "extra  "  rows  from  wrapping  scheme  above 
Fa  =  F(l:cols-l,:); 
for  o=l:Mp-2 

Fa  =  vertcat(Fa,  F(o*cols+l:(o+l)*cols-l,:)); 

end 


%  Create  Vertices  matrix 
a  =  [colx, coly, colz]; 
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%  Create  Figure 

figure  1  =  figure('Color',[l  1  1]); 


%  Create  axes 

axes('Visible',’off, ’Parent’, figurel,’CLim’,[-4e+006  4e+006]); 
view([-83  2]); 

colorbar('FontWeight’,’bold') 

title(plottitle, 'Visible', ’On’, ’FontSize',  14, ’FontWeight’, ’bold') 

%  Call  patch  function 

patch('Vertices',a,'Faces',Fa,'FaceVertexCData',colS,'FaceColor’,’interp’,'FaceLighting’,’goura 

ud') 

%  Uncomment  line  below  and  adjust  numbers  to  set  colorbar  scale 

%  axis('CLim',[-2e+008  2e+008]) 


end 
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Usage  Script 

Load  Variables 

D  =  pt.geometry.D;  %Propeller  Diameter 

N  =  pt.geometry.N;  %Propeller  rotation  rate  [RPM] 

RC  =  pt.design.RC;  %Control  Point  (Section)  Radius  [] 

DR  =  pt.design.DR;  %Control  Point  Radius  Difference  [] 

VSTAR  =  pt.design.VSTAR;  %Total  inflow  velocity  [] 

TANBC  =  pt.design.TANBC;  %Tangent  of  inflow  angle  for  each  Control  Point 

BetalC  =  pt.design.BetalC;  %Ideal  inflow  angle  for  each  Control  Point 

CL  =  pt.design.CL;  %On-design  Lift  Coefficient  for  each  Control  Point 

CD  =  pt.design.CD;  %On-design  Drag  Coefficient  for  each  Control  Point 

CoD  =  pt.design.CoD;  %Chord  Length/Diameter  for  each  Control  Point 

xu  =  pt.geometry.xu;  %Upper  section  x  points  -  Leading  edge  to  trailing  edge 

yu  =  pt.geometry.yu;  %Upper  section  y  points  -  Leading  edge  to  trailing  edge 

xl  =  pt.geometry.xl;  %Lower  section  x  points  -  Leading  edge  to  trailing  edge 

yl  =  pt.geometry.yl;  %Lower  section  y  points  -  Leading  edge  to  trailing  edge 

X3D  =  pt.geometry.X3D;  %x  points  which  create  blade  surface 

Y3D  =  pt.geometry.Y3D;  %y  points  which  create  blade  surface 

Z3D  =  pt.geometry.Z3D;  %z  points  which  create  blade  surface 

UASTAR  =  pt.design.UASTAR;  %Axial  blade  influence  velocity  [] 

UTSTAR  =  pt.design.UTSTAR;  %Tangential  blade  influence  velocity  [] 

Vs  =  pt.input.Vs;  %Design  ship  speed  [m/s] 

Js  =  pt.input.Js;  %Design  Advance  Coefficient 

alpha  =  pt.design.alpha;  %[deg]  Section  angle  of  attack 

theta  =  pi/ 180  *  pt.geometry. theta;  %Pitch  angle  in  radians 

Section  Centroid  and  Moments  of  Inertia  Calculation 

[Mpl,  Ixc,  lye,  Ixyc,  A,  Xbar,  Ybar,  xl,  yl,  xu,  yu]  =  Momenta flnertia(xl,xu,yl,yu); 

Centrifugal  Force 

[omega,  Fc,  Rdif]  =  CentrifugalForce(Mpl,  N,  R,  A,  RC,  DR); 

Stress  Calculation 

[s]  =  Stress(CD,  CL,  BetalC,  Mpl,  Np,  xu,  yu,  xl,  yl,  rho,  Vs,  R,  VSTAR,  CoD,  Rdif,  DR, 
theta,  Xbar,  Ybar,  Ixc,  lye,  Ixyc,  Fc,  A); 

Make  Stress  Plot 

Plot_Blade_Contours(X3D,  Y3D,  Z3D,  s,  'Suction  Side’) 
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Blade  Thickness  ABS 

function  [stress]  =  StructureABS(pt) 

Unpack  Variables 

Datestring  =  pt.date; 
filename  =  pt.  input,  filename; 

Make2Dplot_flag  =  pt.input.Make2Dplot_flag; 

Make3Dplot_flag  =  pt.input.Make3Dplot_flag; 

MakeRhinoflag  =  pt.input.MakeRhinoflag; 

Meanline  =  pt.  input.  Meanline; 

Thickness  =  pt.input.Thickness; 

XR  =  pt.input.XR; 

fOocO  =  pt.input.fOocO; 

tOocO  =  pt.input.tOocO; 

XCoD  =  pt.input.XCoD; 

skewO  =  pt.input.skewO;  %  [deg] 

rakeO  =  pt.input.rakeO; 

alphal  =  pt.input.alphal;  %  [deg] 

CLI  =  pt.  input.  CLI; 

Z  =  pt.input.Z; 

Js  =  pt. input.  Js; 

Vs  =  pt.input.Vs;  %  [m/s] 

R  =  pt.input.R;  %  [m] 

Rhub  =  pt.input.Rhub;  %  [m]  hub  radius 

Np  =  pt.input.Np; 

Hflag  =  pt.input.Hflag; 

Dflag  =  pt.input.Dflag; 

RC  =  pt.design.RC; 

RV  =  pt.design.RV; 

CL  =  pt.design.CL; 

Beta_c  =  atand(pt.design.TANBC);  %  [deg] 

Betal_c  =  pt.design.BetalC*  180/pi;  %  [deg] 

CoD  =  pt.design.CoD; 

TAU  =  pt.design.TAU; 

TANBIV  =  pt.design.TANBIV; 

D  =  2*R;  %  [m] 

Dhub  =  2* Rhub;  %  [m] 

RhuboR  =  Rhub/R; 

N  =  60*Vs/(Js*D);  %  [RPM] 

Interpolate  input  geometry  at  selected  radial  sections 

%  Set  vortex  points  1  and  Mp+1  to  the  hub  and  tip:  RG  =  [0.9*Rhub_oR,RV(2:cnd- 1 ),  1  ]; 

Interpolate  Input  Geometry  at  Sections  with  Cosine  Spacing  Along  the  Span 

RG  =  [0.25,0.70];  %  Required  sections  for  ABS  analysis 
Mp  =  length(RG)- 1 ; 


64 


fOoc  =  pchip(XR,fOocO,RG).*pchip(RC,CL,RG)/CLI;  %  [  ],  camber  ratio  scaled  with  lift 
coefficient 

skew  =  pchip(XR,skewO,RG);  %  [deg],  angular  translation  along  mid-chord  helix 
rake  =  pchip(XR,rakeO,RG)*D;  %  [m],  translation  along  propeller  axis  (3D  X-axis) 

%  tOoc  =  pchip(XR,tOocO,RG);  %  [  ],  thickness  ratio 
%  CoD  =  interpl(RC,CoD,RG,’pchip', 'extrap'); 

CoD  =  pchip(XR,XCoD,RG); 

%  Thesis4:  thickness  profile 

TTRF  =0.5;  %  Tip  Thickness  Reduction  Factor  ==  modified  thickness  at  tip  /  baseline 

thickness  at  tip 

XRmax  =  0.8;  %  maximum  XR  for  which  thickness  reduction  is  less  than  1% 

HTTR  =3.5;  %  Hub-Tip  Thickness  Ratio  =  tO(hub)  /  tO(tip) 

tOtip  =  0.00254;  %  [m]  ==  0.254  mm  =  0.1  inch,  max  thickness  at  tip  section 

tOoc  =  tOtip* (HTTR  -  (HTTR-l).*(RG-Dhub/D)/(l-Dhub/D))./(CoD*D)  .*  (l-(l-TTRF)*exp(- 

4.6*(1-RG)/(1  -XRmax))); 

Find  Basic  Geometry  Parameters  Chord,  Radius,  Pitch,  etc. 

theta  =  interpl(RC,BetaI_c+alphaI.*CL/CLI,RG,'pchip’, 'extrap');  %  Nose-tail  pitch  angle, 
[deg] 

PoD  =  tand(theta).*pi.*RG;  %  Pitch  /  propeller  diameter,  [  ] 
c  =  CoD.*D;  %  section  chord  at  the  sections  [m] 

r  =  RG.*R;  %  radius  of  the  sections  [m] 

theta_Z  =  0:360/Z:360;  %  angle  between  blades  [deg] 

Lay  Out  the  2D  Coordinate  System 

xN  [  ],  x/c  coordinate  in  2D  NACA  foil  tables 
At  the  Leading  Edge:  xN  =  0,  xl  =  c/2,  xO  =  0 
At  the  Trailing  Edge:  xN  =  1,  xl  =  -c/2,  xO  =  1 

xO  [  ],  x/c  distance  along  mid-chord  line  to  interpolate  NACA  foil  table  data. 

xl  [m],  x  distance  along  mid-chord  line  to  evaluate  elliptical  or  parabolic  fonnulae.  By 

definition,  xl  ==  c/2  -  c*x0. 

x2D  [m],  x  position  in  2D  space  on  upper  (x2D_u)  and  lower  (x2D_l)  foil  surfaces 
y2D  [m],  y  position  in  2D  space  on  upper  (x2D_u)  and  lower  (x2D_l)  foil  surfaces 
x2Dr  [m],  x  position  in  2D  space  after  rotation  for  pitch  angle 
y2Dr  [m],  y  position  in  2D  space  after  rotation  for  pitch  angle 

xN  =  [0  .5  .75  1.25  2.5  5  7.5  10  15  20  25  30  35  40  45  50  55  60  65  70  75  80  85  90  95  100]./100; 

xO  =  zeros(l,Np); 
xl  =  zeros(Mp+l,Np); 
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%  %  Even  spacing  along  the  chord 

%  %  for  i  =  1  :Mp  %  for  each  radial  section  along  the  span 

%  for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

%  forj  =  l:Np  %  for  each  point  along  the  chord 

%  xO(lj)  =  (j-l)/(Np-l);  %  [0  :  1] 

%  xl(ij)  =  c(i)/2  -  c(i)*(j-  l  )/(Np-l );  %  [c/2  :  -c/2] 

%  end 
%  end 


%  Cosine  spacing  along  the  chord 

for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

for  j  =  1  :Np  %  for  each  point  along  the  chord 

xl(i,j)  =  c(i)/2  -  0.5*c(i)*(l-cos(pi*(j-l)/(Np-l)));  %  [c/2  :  -c/2] 
end 
end 

xO  =  0.5-xl  ( 1  ,:)/c(l); 

Find  Meanline  and  Thickness  Profiles  (at  xl  positions) 

foe  =  camber  /  chord  ratio  (NACA  data  at  xN  positions) 

dfdxN  =  slope  of  camber  line  (NACA  data  at  xN  positions) 

fscale  =  scale  to  set  max  camber  ratio  to  fOoc  for  each  section 

tscale  =  scale  to  set  max  thickness  ratio  to  tOoc  for  each  section 

f  =  camber  at  xl  positions 

dfdx  =  slope  of  camber  line  at  xl  positions 

t  =  thickness  at  xl  positions 

t  =  zeros(Mp+l,Np); 
f  =  zeros(Mp+l,Np); 
dfdx  =  zeros(Mp+l,Np); 

if  Meanline==0  |  strcmp(Meanline,’NACA  a=0.8  (modified)')  %  Use  NACA  a=0.8  (modified) 
meanline 

foe  =  [0  0.281  0.396  0.603  1.055  1.803  2.432  2.981  3.903  4.651  5.257  5.742  6.120  6.394 
6.571  6.651  6.631  6.508  6.274  5.913  5.401  4.673  3.607  2.452  1.226  0  ]./100; 

dfdxN  =  [0  0.47539  0.44004  0.39531  0.33404  0.27149  0.23378  0.20618  0.16546  0.13452 
0.10873  0.08595  0.06498  0.04507  0.02559  0.00607-0.01404  -0.03537-0.05887-0.08610 
-0.12058  -0.18034  -0.23430  -0.24521  -0.24521  -0.24521]; 

fscale  =  fOoc  /  max(foc); 

dfdxLE  =  0.47539*fscale;  %  slope  at  leading  edge 

for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

for  j  =  l:Np 

f(i,:)  =  pchip(xN,foc  .*fscale(i).*c(i),x0); 
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dfdx(i,:)  =  pchip(xN, dfdxN.  *fscale(i)  ,xO); 

end 
end 

elseif  Meanline==l  |  strcmp(Meanline,'NACA  a=0.8’)  %  Use  NACA  a=0.8  meanline 

foe  =  [0  .287  .404  .616  1.077  1.841  2.483  3.043  3.985  4.748  5.367  5.863  6.248  6.528  6.709 
6.790  6.770  6.644  6.405  6.037  5.514  4.771  3.683  2.435  1.163  0]./100; 

dfdxN  =  [0  .48535  .44925  .40359  .34104  .27718  .23868  .21050 
16892  .13734.11101  .08775  .06634  .04601  .02613  00620  -.01433  -.03611  -.06010  -.08790 - 
.12311  -.18412  -.23921  -.25583  -.24904 -.20385]; 

fscale  =  fOoc  /  max(foc); 

dfdxLE  =  0.48535*fscale;  %  slope  at  leading  edge 

for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

for  j  =  l:Np 

f(i,:)  =  pchip(xN,foc  .*fscale(i).*c(i),x0); 

dfdx(i,:)  =  pchip(xN, dfdxN.  *fscale(i)  ,x0); 

end 
end 

%  elseif  Meanline==2  |  strcmp(Meanline, ’parabolic’)  %  Use  parabolic  meanline 


% 

% 

%  for  i  =  1  :Mp 

%  for  each  radial  section  along  the  span 

% 

for  i  =  1  :Mp+l 

%  for  each  radial  section  along  the  span 

% 

for  j  =  l:Np 

% 

f(i,j)  = 

f0oc(i)*c(i)*(l-(2*xl(i,j)/c(i))A2); 

% 

dfdx(i,j)  = 

-8*f0oc(i)*xl(i,j)/c(i); 

% 

end 

% 

end 

end 

if  Thickness==  1  |  strcmp(Thickness,’NACA  65 A0 10’)  %  Use  NACA  65 A0 10  thickness  form 
toc_65  =  [0  .765  .928  1.183  1.623  2.182  2.65  3.04  3.658  4.127  4.483  4.742  4.912  4.995  4.983 
4.863  4.632  4.304  3.899  3.432  2.912  2.352  1.771  1.188  .604  .021]./100; 

tscale  =  tOoc  /  max(toc_65); 

rLE  =  0.00639*c.*tscale;  %  leading  edge  radius 

for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 
for  j  =  l:Np 

t(i,:)  =  pchip(xN,toc_65.*tscale(i).*c(i),xO); 
end 
end 
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elseif  Thickness==2  |  strcmp(Thickness, 'elliptical')  %  Use  elliptical  thickness  form 
for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

for  j  =  l:Np 

t(i,j)  =  t0oc(i)*c(i)*real(sqrt(l-(2*xl(i,j)/c(i))A2)); 
end 
end 

rLE  =  0;  %  leading  edge  radius 

elseif  Thickness==3  |  strcmp(Thickness, 'parabolic')  %  Use  parabolic  thickness  form 
for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

for  j  =  l:Np 

t(i,j)  =  t0oc(i)*c(i)*(l-(2*xl(i,j)/c(i))A2); 
end 
end 


rLE  =  0;  %  leading  edge  radius 


elseif  Thickness==4  |  strcmp(Thickness,’NACA  65 AO  10  (modified)')  %  Use  modified  NACA 
65 AO  10  thickness  fonn 


xx65mod  =  [0  0 
0.025000000000000 
0.150000000000000 
0.350000000000000 
0.575916230366492 
0.785340314136126 
0.968586387434555 
0.997382198952880 


.005000000000000  0.007500000000000  0. 


0.050000000000000  0 
0.200000000000000  0 
0.400000000000000  0 
0.628272251308901  0 
0.837696335078534  0 
0.981675392670157  0 
1.000000000000000]; 


075000000000000 

250000000000000 

471204188481675 

680628272251309 

.890052356020942 

989528795811518 


012500000000000 

0.100000000000000 

0.300000000000000 

0.523560209424084 

0.732984293193717 

0.942408376963351 

0.994764397905759 


tt65mod  =  [0  0.007650000000000  0.009280000000000  0.011830000000000 


0.016230000000000 

0.036580000000000 

0.049120000000000 

0.046320000000000 

0.029120000000000 

0.008960000000000 

0.004049015364794 


0.021820000000000 

0.041270000000000 

0.049950000000000 

0.043040000000000 

0.023520000000000 

0.007499530848329 

0.000210000000000] 


0.026500000000000 

0.044830000000000 

0.049830000000000 

0.038990000000000 

0.017710000000000 

0.006623639691517 


0.030400000000000 

0.047420000000000 

0.048630000000000 

0.034320000000000 

0.011880000000000 

0.006040000000000 


tscale  =  tOoc  /  max(tt65mod); 


rLE  =  0.00639*c.*tscale;  %  leading  edge  radius 

for  i  =  1  :Mp+l  %  for  each  radial  section  along  the  span 

for  j  =  l:Np 

t(i,:)  =  pchip(xx65mod,tt65mod.*tscale(i).*c(i),x0); 


68 


end 

end 

end 

Find  2D  Unrotated  Section  Profiles 

x2D  [m],  x  position  in  2D  space  on  upper  (x2D_u)  and  lower  (x2D_l)  foil  surfaces 
y2D  [m],  y  position  in  2D  space  on  upper  (x2D_u)  and  lower  (x2D_l)  foil  surfaces 

x2D_u  =  zeros(Mp+l,Np);  x2DI  =  zeros(Mp+l,Np); 
y2D_u  =  zeros(Mp+l,Np);  y2DI  =  zeros(Mp+l,Np); 

for  i  =  1  :Mp+l  %  for  each  section  along  the  span 

for  j  =  1  :Np  %  for  each  point  along  the  chord 

x2D_u(i,j)  =  xl(i,j)  +  (t(i,j)/2)*sin(atan(dfdx(i,j)));  %  2D  upper  surface  x 
x2D_l(i,j)  =  xl(i,j)  -  (t(i,j)/2)*sin(atan(dfdx(i,j)));  %  2D  lower  surface  x 
y2D_u(i,j)  =  f(i,j)  +  (t(i,j)/2)*cos(atan(dfdx(i,j)));  %  2D  upper  surface  y 
y2D_l(i,j)  =  f(i,j)  -  (t(i,j)/2)*cos(atan(dfdx(i,j)));  %  2D  lower  surface  y 
end 
end 

%  Put  all  the  numbers  in  one  list 

%%j  =  l  ==  tail 
%  %  j  =  1  :Np  ==  suction  side 
%  %  j  =  Np  ==  nose 
%  %  j  =  Np  +  1  ==  nose 

%  %  j  =  Np+  1 :2*Np  ==  pressure  side 
%  %  j  =  2*Np  ==  tail 

%  %  Tail  ->  suctioin  side  ->  nose,  nose  ->  pressure  side  ->  tail 

x2D(:,  1  :Np  )  =  x2D_u(:,Np:- 1 : 1);  %  The  first  Np  values  are  the  upper  surface  (suction  side), 

x2D(:,Np+l:Np+Np)  =  x2D_l(:,l:Np);  %  and  the  second  Np  values  are  the  lower  surface 

(pressure  side). 

y2D(:,  l:Np  )  =  y2D_u(:,Np:-l:l); 
y2D(:  ,Np+ 1  :Np+Np)  =  y2D_l(:,l:Np); 

Find  power  at  design  power: 

Pdesign  =  pt.design.CP  *  (1/2)  *  pt.input.rho  *  VsA3  *  pi*RA2; 

Find  expanded  area  ratio 

EAR  =  Z  *  trapz(XR*R,XCoD*D)  /  (pi*RA2); 

Variable  Definitions 

%  a  -  Expanded  blade  area  divided  by  disc  area 
%  a_s  -  Area  of  expanded  cylindrical  section  at  0.25  radius  (mmA2) 

%  C_n  -  Section  modulus  coefficient  at  the  0.25  radius.  To  be  determined 
%  using  the  equation  listed  in  the  code  below. 
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%  C_s  -  Section  area  coefficient  at  0.25  radius  and  is  to  determined  by  the 
%  equation  listed  in  the  code  below. 

%  D  -  Propeller  Diameter  (m) 

%  f,w  -  Material  constants 
%  H  -  Power  at  Rated  Speed  (kW) 

%  I_o  -  Moment  of  inertia  of  expanded  cylindrical  section  at  0.25  radius 
%  about  a  line  through  the  center  of  gravity  parallel  to  the  pitch  line  or 
%  to  the  nose  (mmA4) 

%  K  -  Rake  of  propeller  blade  in  mm  (positive  for  aft  rake  and  negative  for 
%  forward  rake 

%  K1  -  Coefficient  337  for  SI  units 
%  N  -  Number  of  blades 

%  P25R  -  Pitch  at  1/4  radius  divided  by  propeller  diameter, 

%  corresponding  to  the  design  ahead  condition 
%  P25  -  Pitch  at  1/4  radius 

%  P70R  -  Pitch  at  7/10  radius  divided  by  propeller  diameter, 

%  corresponding  to  the  design  ahead  condition 
%  P70  -  Pitch  at  7/10  radius 
%  R  -  rpm  at  rated  speed 
%  S  -  Factor 

%  t25  -  minimum  required  thickness  at  the  thickest  part  of  the  blade 
%  section  at  one  quarter  radius  (mm) 

%  T  -  Maximum  designed  thickness  of  blade  section  at  0.25  radius  from 
%  propeller  drawing  (mm) 

%  theta  -  Pitch  angle  in  radians 

%  U_f  -  Maximum  nominal  distance  from  the  moment  of  inertia  axis  to  points 
%  of  the  face  boundary  (tension  side)  of  the  section  (mm) 

%  W  -  Expanded  width  of  a  cylindrical  section  at  0.25  radius  (mm) 

Material  Data 

%  matl=input('Enter  a  material  type  from  the  list  belowAn  Manganese  bronze\n  Nickel- 
manganese  bronze\n  Nickel-aluminum  bronze\n  Manganese-nickel-aluminum  bronze\n  Stainless 
steel\n  \n','s'); 

%  if  strcmp(matl, ’Manganese  bronze') 

%  f=2.10;  w=8.3; 

%  elseif  strcmp(matl, 'Nickel-manganese  bronze') 

%  f=2.13;  w=8.0; 

%  elseif  strcmp(matl, 'Nickel-aluminum  bronze') 

%  f=2.62;  w=7.5; 

%  elseif  strcmp(matl, ’Manganese-nickel-aluminum  bronze') 

%  f=2.37;  w=7.5;  %  elseif  strcmp(matl, 'Stainless  steel') 

%  f=2.10;  w=7.75;  %  end 

f  =  2.62; 
w  =  7.5; 
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Data  Imported  from  other  OpenProp  Modules 

%  x2D(i,j)  ==  [m]  jth  point  at  ith  blade  section 
%  y2D(i,j)  ==  [m]  jth  point  at  ith  blade  section 
%  i  ==  1  for  r/R  ==  0.25 
%  i  ==  2  for  r/R  ==  0.70 

%  Tail  ->  suctioin  side  ->  nose,  nose  ->  pressure  side  ->  tail 

%  j  =  1  ==  tail 

%  j  =  1  :Np  ==  suction  side 

%  j  =  Np  ==  nose 

%  j  =  Np  +  1  ==  nose 

%  j  =  Np+  1 :2*Np  ==  pressure  side 

%  j  =  2*Np  ==  tail 

%  Pdesign  ==  [W]  power  required  at  the  design  state 
%  Z  ==  number  of  blades 

%  D  ==  [m],  diameter 

%  N  ==  [RPM],  rotation  rate 
%  EAR  ==  [  ],  expanded  blade  area  /  disc  area 

%Points  of  foil  section  at  r=0.25R  and  r=0.70R 
for  n=l:2 
if  n==  1 

xu25  =  x2D(n,l:Np)*1000; 
xl25  =  x2D(n,Np+l:2*Np)*1000; 
yu25  =  y2D(n,l:Np)*1000; 
yl25  =  y2D(n,Np+l:2*Np)*1000; 


elseif  n==2 

xu70  =  x2D(n,l:Np)*1000; 
xl70  =  x2D(n,Np+l  :2*Np)*  1000; 
yu70  =  y2D(n,l:Np)*1000; 
yl70  =  y2D(n,Np+l  :2*Np)*  1000; 

end 

end 

%Flip  lower  surface  points  for  use  in  determining  maximum  blade  thickness 

xl25  =  fliplr(xl25); 
yl25  =  fliplr(yl25); 
xl70  =  fliplr(xl70); 
yl70  =  fliplr(yl70); 

%Power  at  rated  speed  (kW) 

H  =  Pdesign/1000; 

%RPM  at  rated  speed 

R  =  N; 
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%  Number  of  blades 

N  =  Z; 


%Area  ratio 

a  =  EAR; 


%Rake  (mm) 

K=  0; 

%Pitch  ratio  at  1/4  and  7/10  radius 

P25R=PoD(l); 

P70R=PoD(2); 

%Expanded  width  of  1/4  radius  section  (mm) 

Imax25  =  flnd(xu25==max(xu25)); 

Imin25  =  find(xu25==min(xu25)); 

W  =  sqrt((xu25(Imax25)-xu25(Imin25))A2  +  (yu25(Imax25)-yu25(Imin25))A2); 
%  a_s  (mmA2) 

%Shift  of  section  to  accomodate  lower  section  zero  crossings 

ymin25  =  min(yl25); 
yu25  =  yu25-ymin25; 
yl25  =  yl25-ymin25; 

xmin25  =  min(xl25); 
xu25  =  xu25  +  xmin25; 
xl25  =  xl25  +  xmin25; 

Section  Area  and  Centroid  Calculation 

Au  =  trapz(xu25,yu25); 

A1  =  trapz(xl25,yl25); 
a_s  =  (Au  -  Al); 

%  Io  (mmA4) 

%Centroid  Calculation 

hlu=zeros(l,length(yu25)-l); 

hll=hlu; 

for  n=  1  :length(yu25)- 1 
h  Iu(n)=min(yu25(n),yu25(n+ 1)); 
end 

for  n=l:length(yl25)-l 
hll(n)=min(yl25(n),yl25(n+l)); 
end 
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h2u=abs(diff(yu25)); 

h21=abs(diff(yl25)); 

dxu25=diff(xu25); 

dxl25=diff(xl25); 

Mxu25=sum(dxu25/2.*(h2u.A2/3  +  hlu.*h2u  +  hlu.A2)); 
Mxl25=sum(dxl25/2.*(h21.A2/3  +  hll.*h21  +  hll  A2)); 

Mx=Mxu2  5  -Mxl2  5 ; 
ybar=Mx/a_s; 

Moment  of  Inertia  Calculation 

Iu=sum(dxu25.*(4*hlu.A3  +  h2u.A3  +  6*hlu  A2.*h2u)); 

Il=sum(dxl25.*(4*hll.A3  +  h21.A3  +  6*hll  A2.*h21)); 

Ix=Iu-Il; 

Io=(Ix-yb  arA2  *  a_s); 

plot(xu25,yu25,’-bs',xl25,yl25,'-rsVMarkerSize',8,'LineWidth',2) 
grid  on 

x=[min(xu25),max(xu25)] ; 
y=[ybar,ybar]; 

line(x,y,'ColorVg','LineWidth',2,'LineStyleV— ') 

Maximum  Distance  from  Neutral  Axis  to  Tension  Side  (mm) 

U_f=max(abs(ybar-yl25)); 

Maximum  Designed  Blade  Thickness  (mm) 

no  =  round(Np*.25); 
thick  =  zeros(l,Np-no); 
for  n=Np:-l:no; 

thick(n)=abs(yu25(n)-yl25(n)); 

end 

T  =  max(thick); 

[trash,Ithickmax]  =  max(thick); 
x  =  [xu25(Ithickmax),xu25(Ithickmax)]; 
y  =  [yl25(Ithickmax),yu25(Ithickmax)]; 
line(x,y,'Color','rVLineWidth',2,'LineStyleV— ') 
xlim(  [min(xu25)- 1 0,max(xu25)+ 1 0]); 
ylim(  [min(yl25)- 1  ,max(yu25)+ 1  ]); 

%  Iotext=num2str(Io); 

text(xmin25,T/2,strcat(’Io  =  ’,num2str(fix(Io)),’  mmA2’),'fontweighf,’b’,'fontsize',12) 
%  annotation('textbox', [0 .25 ,0 .25 , .3 ,.3] ) 
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title(’Section  at  25%  RVfontsize',14,'fontweight’,'b') 
xlabel(’X2D  (mm)7fontweight','b') 
ylabel('Y2D  (mm)7fontweight','b') 

legend('suction  side’,'pressure  side',  'neutral  axis',  'max  thickness') 

Determining  S  Factor  for  SI  and  MKS  units 

if  D<=6.1 
S=l; 
else 

S=sqrt((D+24)/30. 1); 
end 

if  S>  1.025 
S=1.025; 
end 

Calculation  of  Constants  and  Required  Blade  Thickness 

Kl=337;  %Constant  for  SI  and  MKS  units 
A=  1 +6/P70R+4 . 3  *  P2  5  R; 

B=(4300*w*a/N)*(R/100)A2*(D/20)A3; 

C=(l+1.5*P25R)*(W*f-B); 

C_n=Io/(U_f*W*TA2); 

C_s=a_s/(W*T); 

t25=S*(Kl*sqrt(A*H/(C_n*C*R*N))+(C_s/C_n)*B*K/4/C); 

Output  stress  data  structure 

stress. Pdesign  =  Pdesign; 
stress. Z  =  N; 
stress. D  =  D; 
stress  .N  =  R; 
stress. EAR  =  EAR; 
stress.  t25  =  t25; 
stress. T  =  T; 
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Automation  Direct _ Braking  Resistor,  1000  W,  50  Ohm,  NEMA  1  Enel. _ GS-4015-BR-ENC _ 1 _ $219.00 _ $219.00 

FilterConcepts _ EMI  Line  Filter _ SF20L _ 1 _ S86.13 _ S86.13 

McMaster  Carr _ Long  Fastener _ Transition  Motor  Housing  Fastener _ 92196A287 _ 3 _ $3.68 _ $11.04 

McMaster  Carr _ Customized  -  Turned  Head _ Brush  Plate  Fastener _ 93996A531 _ 2 _ $3.65 _ $3.65 

McMaster  Carr _ Insert  for  nose  fairing _ Threaded  Insert _ 92066A036 _ 1 _ S4.92 _ S4.92 
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